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SUMMARY 


The  development  of  a tool  for  solving  the  near  field  of  a scram- 
jet  fuel  injector  was  attacked  by  first  developing  a numerical  technique 
for  solving  the  laminar,  supersonic  near  wake  flow.  It  was  considered 
important  to  develop  a procedure  that  had  a potential  for  reduced  compu- 
tation time  compared  with  explicit  methods.  The  implicit  numerical 
procedure  «f-Rr±iey-ana  McDonald  was  extended  to  mixed  subsonic/super- 
sonic flow  with  shocks,  expansions,  and  regions  of  reverse  flow.;  Briley 
and  McDonald  had  previously  applied  the  procedure  to  subsonic,  constant 
area  duct  flow  with  no  recirculation.^  In  the  present  case,  numerical 
results  have  been  obtained  for  the  laminar,  supersonic  near  wake  behind 

-4 

a rectangular  base.  f~*  \ * /.V 

The  numerical  method  applies  a time  linearization  based  on  a 
Taylor  series  expansion  about  the  known  time  level,  and  the  Douglas-Gunn 
alternating  Direction  Implicit  (ADI)  procedure  to  the  Navier-Stokes 
equations.  Briley  and  McDonald  obtained  the  finite  difference  equations 
by  using  standard  three-point  central  differencing.  This  generated  a 
series  of  block  tridiagonal  systems  which  can  be  quickly  solved  by  a 
standard  elimination  technique.  The  same  approach  was  followed  here, 
except  that  all  of  the  differential  equations  were  written  in  the  conser- 
vation form  (Briley  and  McDonald  used  the  non-conservation  form  of  the 
energy  equation)  and  the  finite  difference  equations  were  derived  by  the 
cell  integration  technique.  The  cell  Integration  technique  considers 
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the  conservation  equations  as  integral  laws  over  a control  volume  (cell) 
around  a grid  point  and  also  leads  to  central  differencing  for  the 
interior  grid  points. 

The  chief  advantage  in  using  the  cell  integration  technique  is 
the  conceptual  aid  afforded  in  applying  the  boundary  conditions.  Allen 
and  Cheng  used  this  technique,  and  their  work  served  as  a guide  in 
selecting  one-sided  difference  forms  for  the  nonzero  boundary  terms. 
Because  the  present  method  is  implicit,  however,  several  new  forms  were 
required  for  stable  and  accurate  solutions.  It  was  found,  for  example, 
that  second-order  forms  for  the  pressure  and  3v/3y  on  the  centerline  are 
needed  to  prevent  y-direction  wiggles  in  the  steady  state  solution.  Also, 
a new  implicit,  linear  extrapolation  scheme  using  the  finite  difference 
equations  was  developed  for  the  outflow  boundary.  This  was  required  to 
eliminate  wiggles  in  the  x-direction  in  the  steady-state  solution.  All 
the  explicit  extrapolation  schemes  at  the  outflow  caused  the  solution 
to  diverge  for  At  > At„_ . Zero-gradient  forms  at  the  outflow  boundary, 

Lr  L 

whether  explicit  or  implicit,  caused  x-direction  wiggles  in  the  steady 
state  solution. 

Three-dimensional  contour  plots  proved  to  be  an  important  diag- 
nostic tool.  It  was  not  discovered  that  the  x-direction  wiggles  were 
caused  by  the  treatment  of  the  outflow  boundary  conditions  until  the 
3-D  plots  clearly  revealed  that  as  the  recompression  wave  crossed  the 
downstream  boundary,  the  wiggles  formed  and  propagated  upstream  to  the 
back  wall  and  inflow  regions.  Up  to  then,  the  suspected  causes  were 
improper  treatment  of  the  back  wall  boundary  conditions,  or  that  the 
cell  Reynolds  numbers  were  greater  than  two  there.  The  use  of  upwind 
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differencing,  artificial  viscosity,  or  a much  smaller  Ax  were  considered 
to  be  undesirable  remedies. 

. The  results  for  the  contour  plots  showed  qualitative  agreement 

with  Allen  and  Cheng  and  Kronzon,  et  al.,  and  close  quantitative  agree- 
ment where  comparisons  were  possible.  The  centerline  pressure  plot 
showed  very  close  quantitative  agreement  with  Allen  and  Cheng.  As  a 
further  check  on  accuracy,  overall  mass  balances  were  computed  at  each 
time  step.  In  the  (nearly)  steady  state  conditions,  net  mass  inflow 
rate  differed  from  net  mass  outflow  rate  by  about  1.8%  or  less. 

No  artificial  viscosity  was  used  in  obtaining  these  solutions. 

It  is  interesting  that  Briley  and  McDonald  required  an  additional 
explicit  artificial  viscosity  in  their  subsonic  duct  flow  solutions. 

The  reasons  for  this  difference  in  behavior  are  not  known.  It  may  be 
speculated,  however,  that  the  difference  arises  from  the  present  use  of 
the  conservative  form  of  the  conservation  equations,  the  cell  integration 
technique  for  generating  finite-difference  equations,  and  the  correspond- 
ing careful  treatment  of  the  boundary  conditions. 

A time  step  limitation  was  expected,  although  the  method  is 
• implicit,  because  the  equations  were  linearized  with  respect  to  time. 

For  one  set  of  initial  conditions,  this  limitation  was  found  to  be  around 
32  Atgp^.  The  present  method  had  a computation  time  per  time  step  per 
grid  point  of  approximately  five  times  longer  than  Allen's  explicit 
method,  but  could  take  time  steps  over  30  times  larger.  This  represents 
a six-fold  decrease  in  computation  time.  In  addition,  the  ability  to 
change  (increase)  the  size  of  the  time  step  during  computation  to  reduce 
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computation  time  was  demonstrated.  This  suggests  that  a time  step 
strategy  might  be  successful  wherein  smaller  At's  were  used  at  the 
beginning,  followed  by  increasing  At  as  the  steady  state  is  approached. 
This  would  be  appropriate  when  the  assumed  initial  conditions  were  very 
far  from  the  steady-state  solution.  Thus  the  method  appears  to  offer 
significant  time  savings. 

The  effect  of  initial  conditions  on  the  steady-state  solution 
was  examined.  To  do  this,  a range  of  initial  horizontal  velocities  were 
applied  in  the  region  below  the  expansion  corner.  All  the  other  initial 
conditions  were  the  same:  a boundary  layer  on  the  upper  wall  upstream 
of  the  expansion  corner  and  freestream  conditions  elsewhere.  It  was 
shown  that  u « 0 led  to  divergence  for  At  * 16AtCF^.  Increasingly 
rapid  rates  of  convergence  were  realized  as  u was  increased  from  10%  to 
100%  of  the  freestream  value.  The  results  for  all  the  converged  cases 
indicated  that  the  final  solution  was  insensitive  to  the  initial  condi- 
tions, but  that  the  time  to  convergence  was  highly  dependent  on  initial 
conditions.  Also,  convergence  was  shown  to  occur  for  a significant 
range  of  initial  backwall  u. 

Accuracy  of  the  coarse  mesh  results  was  shown  by  comparisons  with 
the  fine  mesh  solution.  Both  solutions  were  in  close  agreement.  Small, 
irregular  disturbances  in  the  inflow  region  and  in  the  shock  occurred 
for  the  coarse  mesh  solution.  These  can  be  attributed  to  the  lack  of 
resolution  in  the  coarse  mesh  in  the  Inflow  boundary  layer  and  in  the 
shock  at  the  outflow,  as  they  disappeared  in  the  fine  mesh  solution. 
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These  numerical  results  served  to  demonstrate  that  this  numeri- 

* * 

cal  method  produced  stable,  convergent,  and  accurate  solutions  when 
applied  to  this  complex  flow  problem. ^To^the  author's  knowledge,  no 
other  implicit  scheme  has  been  successfully  applied  to  the  multidimen- 
sional nonlinear  Navier-Stokes  equations  for  the  supersonic  base  flow 
problem. 
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CHAPTER  I 

INTRODUCTION 

Recent  interest  in  hypersonic  flight  has  motivated  an  increasing 

number  of  investigations  into  advanced  airbreathing  propulsion  devices, 

including  supersonic-combustion  ramjets  (or  scramjets).  Many  studies 

have  been  related  to  an  airbreathing  launch  vehicle  for  NASA's  space 

shuttles,  but  found  that  the  technological  state-of-the-art  of  the  pro- 
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pulsion  system  was  not  sufficiently  developed.  More  recently  atten- 
tion has  been  given  to  developing  a scramjet  engine  for  a hypersonic 
4 

research  vehicle.  A principle  requirement  of  the  scramjet  is  the  speci- 
fication of  the  flow  field  downstream  of  the  fuel  injector.  Knowledge 
of  the  combustion  flow  field  and  heat  release  distributions,  for  example, 
would  allow  for  the  design  of  engines  requiring  a fraction  of  the  fuel 
heat  sink  capacity  for  cooling.  This  would  allow  the  airframe  designer 
more  flexibility.  Additionally,  there  is  the  need  for  complete  combustion 
in  as  short  a distance  as  possible,  so  that  long  combustors  will  not  be 
required.  Hence  the  need  for  rapid  mixing  of  the  fuel  and  air  streams 
makes  the  near  field  of  the  injector  a region  of  great  interest. 

The  injector  flow  field  is  quite  complex,  which  greatly  hinders 
analysis  (see  Figure  1) . Shocks,  high  transverse  pressure  gradients, 
and  region  of  reverse  flow  make  the  near  field  similar  to  a base  flow, 
but  with  the  added  complications  of  fuel  injection  and  subsequent  mixing 

and  combustion.  All  of  the  flow  features  strongly  influence  the  turbu- 
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lent  mixing  and  combustion.  Hence  an  accurate  analytical  procedure 
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which  describes  the  fuel  injector  near  field  must  include  all  these 
features.  Previous  studies^-1"*  examined  the  mixing  and  combustion  of 
compressible  turbulent  streams  but  neglected,  and  indeed  could  not  com- 
pute, the  important  effects  of  shocks,  recirculation,  and  regions  of 
high  transverse  pressure  gradients. 

The  present  analysis  seeks  to  include  these  effects,  but  neglects 
the  turbulence  and  combustion  for  two  reasons.  Uncertainties  in  the 
turbulence  and  combustion  models  limit  the  validity  of  analysis.  Even 
relatively  advanced  turbulence  models,  such  as  those  where  velocity  and 
length  scales  are  computed  from  differential  equations, ^ ^ may  have 
difficulty  in  describing  details  of  this  flow.  Thus,  it  would  be  diffi- 
cult to  establish  whether  inaccuracies  in  a new  numerical  procedure  were 
due  to  the  models  or  to  the  method  itself.  Second,  experience  has  shown 

the  differential  equations  of  turbulence  to  be  troublesome  numerically, 

20 

which  greatly  hinders  even  the  development  of  a new  method.  It  is 

therefore  prudent  to  prove  a new  method  first  by  solving  a similar 

problem  where  the  flow  is  well-characterized  by  the  governing  equations. 

Here,  the  laminar,  supersonic  base  flow  problem  (Figure  2)  was  solved. 

The  flow  is  specified  by  the  Navier-Stokes  equations  along  with  the 

conservation  equations  of  mass  and  energy  and  the  equation  of  state. 

Many  previous  base  flow  studies  used  an  integral  technique  to 

determine  a base  pressure  or  a base  drag,  but  few  other  details  of  the 
21 

flow.  Mueller,  for  example,  used  the  Chapman-Korst  method  to  determine 
a single  turbulent  base  pressure  for  supersonic  axlsymmetric  flow,  which 
was  assumed  to  be  constant  across  the  base.  Mueller  pointed  out  that 
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Figure  2.  No  Injection  Flow  Field  (Base  Flow) 
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the  solutions  obtained  were  asymptotic  and  valid  only  for  high  Reynolds 

numbers,  and  that  the  important  effects  of  the  initial  boundary  layer 

22 

were  neglected.  Alber  and  Lees,  used  the  Crocco-Lees  integral  proce- 
dure for  supersonic  turbulent  base  flows.  They  showed  that  the  initial 
boundary  layer  can  dominate  the  viscous  Interaction  near  the  base  when 
its  height  is  of  the  order  of  the  base  height.  They  also  computed  a 
constant  pressure  up  to  the  rear  stagnation  point,  the  distance  to  the 
rear  stagnation  point,  the  centerline  pressure  distribution  downstream 
of  this  point,  and  the  correct  trend  of  increasing  average  base  pressure 
for  increasing  initial  boundary  layer  thickness.  Both  integral  theories 
rely  heavily  on  the  flow  being  well-characterized,  beforehand  by  another 
method  or  by  experimental  data.  Neither  have  been  shown  to  compute 
details  within  the  recirculation  region,  variation  of  pressure  within 
the  recirculation  region  and  along  the  base,  and  shocks.  Extension  of 
an  integral  technique  to  include  these  features  and  subsequent  extension 
to  the  case  of  a fuel  injector  appears  to  be  unpromising  at  best.  Even 

extension  of  the  Crocco-Lees  method,  for  example,  to  axisymmetric  flow 

23 

has  been  accomplished  only  with  great  difficulty  by  Mehta. 

Finite  difference  procedures  appear  to  be  more  promising  for 

computing  the  base  flow  field.  So  far  as  the  author  knows,  the  only 

finite  difference  techniques  applied  to  the  full  conservation  equations 

for  the  supersonic  base  flow  problem  were  the  explicit  methods  of  Allen 
24  25 

and  Cheng  and  of  Roache  and  Mueller.  For  the  latter  case,  however, 
relatively  little  Information  about  the  solution  was  provided.  The 
Allen  and  Cheng  method  computed  the  steady  state  solution  by  solving  the 
unsteady  equations  for  asymptotically  large  time.  However,  much 
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computation  time  was  required  because  the  method,  as  all  explicit  methods, 

is  subject  to  one  or  more  stability  restrictions  on  the  time  step  size 

relative  to  the  spatial  grid  size.  These  stability  criteria  are  the 

well  known  Courant-Friedrichs-Lewy  (CFL)  condition  (in  one  dimension, 

AtCFL  — Ax/(|u|  + c))  and,  in  some  methods,  a viscous  stability  limit 
2 

(At  <_  Ax  /2v).  Since  the  maximum  time  step  size  is  related  to  the  spa- 
tial grid  size,  when  accuracy  is  desired  and  a fine  mesh  is  used,  the 
computation  time  correspondingly  increases. 

Implicit  methods,  on  the  other  hand,  tend  to  be  stable  for  much 
larger  time  steps.  Hence  they  offer  the  prospect  of  faster  solution 
than  explicit  methods,  provided  the  computation  time  per  time  step  is 
comparable  to  that  of  explicit  methods.  When  applied  to  one-dimensional 
equations  using  central  differencing,  an  implicit  method  usually  gives 
a linear  system  with  a tridiagonal  coefficient  matrix  which  is  easily 
and  quickly  solved.  Multidimensional  problems,  however,  give  more  com- 
plicated coefficient  matrices  which  are  time  consuming  to  solve.  Further 

the  equations  need  to  be  suitably  linearized  before  application  of  the 

26 

implicit  technique.  Briley  and  McDonald  have  proposed  a procedure  which 
linearizes  the  unsteady  equations  in  time  by  Taylor  series  expansion  about 
the  known  time  level.  It  preserves  the  efficiency  of  one-dimensional 
systems  by  applying  an  Alternating-Direction- Implicit  (ADI)  procedure, 
in  which  the  equations  are  considered  implicit  in  one  direction  at  a 

27 

time.  The  particular  ADI  scheme  used  here  is  that  of  Douglas  & Gunn. 

This  method  is  tentative  because  Briley  and  McDonald  only  applied 
the  method  to  a subsonic  duct  flow  with  no  recirculation.  This  is  very 
different  from  the  supersonic  base  flow  problem,  and  the  ability  to 
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compute  shocks  and  recirculation,  for  example,  needs  to  be  proven.  In 
addition,  Briley  and  McDonald  never  established  the  accuracy  of  the 
method  by  comparing  with  experimental  data  or  an  exact  solution.  At  most 
they  have  only  shown  qualitative  agreement  with  approximate  (i.e.,  one- 
dimensional exact)  analyses.  Nevertheless,  the  method  appears  to  be 
promising  in  not  being  subject  to  stability  limits  on  the  time  step 
size  and  in  retaining  the  computational  speed  of  one-dimensional  implicit 
systems. 

26 

In  this  thesis,  the  Briley  and  McDonald  procedure  was  applied 
to  the  governing  equations,  and  the  cell  integration  technique  was 
applied  to  derive  the  finite  difference  equations.  In  brief,  the  govern- 
ing equations  were  linearized  in  time  by  a Taylor  series  expansion  about 
the  known  (or  nC^)  time  level.  The  finite  difference  equations  were 
derived  by  applying  the  cell  integration  technique,  which  leads  to 
central  differencing  for  the  interior  grid  points.  Application  of  the 
ADI  procedure  leads  to  sequences  of  one-dimensional  implicit  systems 
having  block  tridiagonal  coefficient  matrices.  Each  of  these  systems 
(one  sequence  of  systems  for  each  coordinate  direction)  is  solved  by 

the  standard  block  elimination  technique  as  outlined  in  Isaacson  and 

28 

Keller.  No  iteration  is  required  to  compute  the  solution  for  a given 
time  step. 

The  method  was  checked  against  the  previous  laminar,  supersonic 

24  29 

base  flow  calculations  of  Allen  and  Cheng  and  Kronzon,  et  al.  This 
allowed  a check  of  the  capability  of  this  implicit  method  to  compute  a 
flow  with  shocks,  reverse  flow,  high  transverse  pressure  gradients,  and 
a wide  variety  of  boundary  conditions.  It  also  allowed  a check  on  the 
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ability  of  the  method  to  compute  a solution  in  less  time  than  by  an 
explicit  method. 
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CHAPTER  II 

GOVERNING  EQUATIONS 

The  governing  equations  are  the  conservation  equations  for  mass, 

x-momentum,  y-momentum,  and  energy,  and  the  equation  of  state  for  the 

two  dimensional  flow  of  a perfect  gas  with  constant  specific  heats. 

The  differential  equations  are  written  in  the  conservation  form.  As 
30 

Roache  shows  (p.  28),  when  the  conservation  form  is  used,  then  the 
finite  difference  equations  preserve  the  integral  Gauss  divergence 
property  of  the  continuum  equation.  His  example  illustrates  the  alge- 
braic balance  of  flux  quantities  and  accumulation  rates  in  a small  con- 
trol volume.  This  has  an  intuitive  appeal.  In  addition,  Roache  points 
out  that  the  Rankine-Hugoniot  relations  were  derived  from  the  conserva- 
tion form  and  hence  the  jump  conditions  across  a shock  are  automatically 
satisfied.  No  special  treatment  is  given  for  the  formation  of  shocks  if 
they  develop.  This  is  called  "shock-capturing"  or  "shock-smearing." 

(In  the  Russian  literature  this  is  called  a "through"  method.  See 
30 

Roache  p.  227). 

The  dimensional  equations  are: 


ft  * ~ lx (pu)  - k (pv) 


k [pu2  + P - Txx]  ' 17  [puv  ' TxyJ 


(2-1) 


3 (pu) 
3t 


(2-2) 
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9(pv)  = _ 1_ 
at  ax 


[puv  - T ] - [pv^  + P - T ] 
3x  k xyJ  3y  * yyJ 


(2-3) 


2 . 2. 


h tpe  + 2 <u2  + v2>^  * ‘ h [pu(e  + p + ~U-p/— ) ~ qj  - |~[pv(e  + 


_£  + _(.M— t V _ n ] + (ui  + VT  ) + (UT  + VT  ) 

p 2 7 3x  xx  xy  3y  xy  yy 


where: 


, 3v.  3v. 

(2-4) 

p - (y  - 1)  e 

(2-5) 

,4  3u  2 3v, 
xx  = u '3  3x  " 3 3y  ; 

(2-6) 

,3u  , 3v. 
xy  ” u ^3y  3x^ 

(2-7) 

,4  3v  2 3u  . 
yy  * M3  3y  ” 3 3x  ; 

(2-8) 

yk  3e 

Sc  c 3x 

(2-9) 

yk  3e 

c 3y 
P 


(2-10) 


Here  c^,  y,  k,  and  p are  assumed  to  be  constant.  (The  case  of  temperature 
variation  of  these  quantities  Is  a straightforward  extension  if  done 
explicitly.)  Bulk  viscosity  was  assumed  to  be  zero. 

Great  utility  is  afforded  by  nondlmensionalizing  the  equations. 


Then,  different  flow  conditions  can  be  characterized  easily  by  a small 
set  of  nondimens lonal  parameters  rather  than  having  to  change  all  the 
dimensional  parameters.  The  following  convention  was  used: 
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where  H ■ base  half  height  and  subscript  1 denotes  freestreara  quanti- 
ties. 

Substituting  these  quantities  into  the  conservation  equations, 
eliminating  pressure  with  the  state  equation,  and  dropping  the  overbars 
equations  (2-1)  through  (2-4)  become 


3 3 / x 3 / \ 

_ . . _ (pu)  . _ (pv) 


(2-11) 


3pu  3 r„  2 / 1 \ 1 ,4  3u  2 3vu  3 , 1 ,3u  3vN1 

-JJ -5J  [pu  + (-j)  pe  - — (T  - - T _ (puv  - jj)] 


(2-12) 


3pv  _ 3 1 x3u  3v..  3 ,2  /lx  1 ,4  3v  2 3u. . 

3t  " 3x  [p  " Re  3y  + 3x)  ^ " 3y  [p  + ( M2)pe  “ Re(3  3y~  3 3x^ 

3 3 


(2-13) 


3_ 

3t 


pK,  2 2. 

~2-(u  + v ) 


3_ 

3x 


3_ 

3y 


t — VUl 

3y  xy 


K . 2 2. 

j (u  +v  ) 


yy 


3u 


where: 


xx 


(u2  + V2))  + yqx] 

[1—  (uT  + VT  ) 

3x  xx  xy 

u . 3Vx  . 3v,  t 

y + IS*  + Tyy  3y]} 

(2-14) 

u 2 3v. 
x " 3 3y; 

(2-15) 

_ - -L  f3u  . 3V) 

Txy  Re  ' 3y  T ax'' 


^ 1 ,4  3v  2_  3u. 

Tyy  “ Re  S 3y  ‘ 3 3x; 


q i — l£ 

nx  RePr  3x 

1_  3e 

^y  RePr  3y 


(2-16) 


(2-17) 


(2-18) 


(2-19) 


K - y(y-1)M.‘ 


plUlH 


1 Y(rl)e, 


These  are  the  equations  solved  in  the  program  along  with  boundary 

* 

conditions  which  are  derived  for  the  geometry  illustrated  in  Figure  2. 

The  finite  difference  grid  for  the  flow  field  is  shown  in  Figure  3. 

Note  that  only  half  of  the  flow  field  is  computed  (the  upper  half)  since 
the  flow  is  assumed  symmetric  about  the  central  plane  (DE) . The  flow 
is  two-dimensional,  planar  over  the  rectangular  corner  BCD.  The  incoming 
flow  is  supersonic  with  a boundary  layer  on  the  upper  wall.  Details  of 
the  boundary  and  initial  conditions  are  given  in  Chapter  5. 


i 


Figure  3.  Flow  Field  Boundary  and  Grid  (not  to  scale). 
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CHAPTER  III 


NUMERICAL  PROCEDURE 


The  unsteady  equations  are  parabolic  in  time  and  at  each  time 
step  the  simultaneous  solution  of  the  equations  for  the  whole  flow 
field  is  required.  Numerically,  this  means  the  solution  of  a linear 
system  (for  implicit  methods)  which,  in  turn,  requires  the  lineariza- 
tion of  the  non-linear  terms  at  the  implicit  time  level.  This  is 
accomplished  by  a Taylor  series  expansion  about  the  solution  at  the 
known  time  level  as  outlined  in  Briley  and  McDonald.  They  point 
out  that  this  procedure,  adapted  for  the  integration  of  initial-value 
problems,  permits  the  computationally  efficient  solution  of  coupled, 
non-linear  equations  in  one  space  dimension  by  a one-step  non-iterative 
scheme.  The  efficiency  is  retained  for  multidimensional  problems  by 
using  alternating-direction  implicit  (ADI)  techniques.  As  an  example  of 
the  linearization,  the  continuity  equation  (2-11)  becomes: 


(Dn+1-Dn'>  3 , ,n+l  3 , vn+1 

at  ? m ‘ 3x  (pu)  * 3^ 


3 rt  \n.  . 3p  . 3u.  3u.n  . 3 r / \n.  3v.n.  . 

“ Jx  + a + + P P At]  - (pv)  + p — ) At] 


3t 


3 1 


3t 


3y 


3t 


or 


, n+1  n. 

ifi ~P  > » 

At 


3 r n+1  n , n n+1  n n,  3 r n+1  n . n n+1  n n. 
- —I  P u+pu  - p u ] - — L p v+pv  -pv] 


3x 


(3-1) 
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Equation  (3-1)  is  linear  in  the  unknown  (or  n+1)  variables.  The  complete 
set  of  linearized  equations  is  given  in  Appendix  A. 

To  obtain  the  finite  difference  form  of  the  equations,  the  cell 
integration  technique  is  used.  This  technique  is  best  Illustrated  by 
example.  Consider  the  differential  equation  for  conservation  of  mass, 
using  Cartesian  tensor  notation: 

3p/3t  + 3pu1/3xi  - 0 (3-2) 

This  equation  is  integrated  over  a control  volume  (cv)  which  is  a cell 
of  dimensions  Ax,  Ay,  Az.  Thus 

///  (3p/3t)dv  = - fff  (3pu  /3x  )dv 
cv  cv  1 1 

■ - //  pu±nidA  (3-3) 

cs 

using  Gauss'  theorem.  For  two-dimensional  flow,  the  area  integral  over 
the  control  surface  becomes: 

//  pu^n^dA  ■ //  pu  dydz  + ff  pu(-l)dydz 
cs  x+  x- 

+ //  pv  dxdz  + //  pv(-l)  dxdz  (3-4) 

y+  y- 

where  x+,  x-,  y+,  y-  are  the  control  surface  or  cell  edges  (see  Figure  4). 
Next,  the  following  approximations  are  made: 
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///  (3p/3t)dv  i (3p/3t)ij  Ax  AyAz  ■ t 

cv 

jj  pu  dydz  < (pu)^  AyAz  (3-5) 

x+ 

etc. 

Equation  (3-4)  then  becomes 

(3p/3t)ij  - - ^ [(P«)x+-  (pu)xJ  ~'^t(pv)y+"  <Pv)yJ 
or  ■ 

3p/3t  «*  -6x(pu)  - <5y(pv)  (3-6) 

when  (pu)x+,  for  example,  is  taken  as  the  average  between  (pu)„  and 

(pu)  _ . and  similarly  for  (pu)  , then  3x(pu)  becomes  the  standard 
i+i , j x— 

central  difference  form.  The  chief  advantage  in  using  this  formulation 
is  the  conceptual  aid  afforded  in  applying  the  boundary  conditions. 

With  this  technique  all  spatial  derivatives  are  related  to  values  at  the 
cell  edges.  Thus,  when  boundaries  are  adjacent  to  cell  edges  it  becomes 
clear  which  terms  must  be  modified  to  match  the  boundary  condition. 

The  complete  set  of  finite  difference  equations  is  given  in  Appendix  B. 

For  interior  points,  the  value  of  a quantity  on  a cell  edge  is  always 
taken  as  the  average  of  the  cell  points  on  either  side  of  the  edge. 

A derivative  at  the  cell  edge  is  always  the  difference  of  the  two  near- 
est cell  points.  Thus,  for  example 

i 

| (pu)x+ " i[^u)ij + (pu)i+i,j]  (3_7) 

lh  (pu)1x+“  zi  [(pu)i+i,j-  (pu)ij] 

t 


(3-8) 
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When  these  are  carried  out  for  all  the  cell  edges,  the  result  is  central 
differencing.  For  example 


. , . 1 Ir(pu)i+ii+(pu)i1  (pu)11+(pu)i-li 

6x(pu)  = — [(pu)^-  (pu)xJ  - — [ 2 2 


2&x  C(pu)i+lj  " (pu)i-lj] 


then  equation  (3-1)  can  be  written  as: 


i>n+1 
^ i+lj 


+ b 


n,n+l 

ij+1 


+ c 


n ,n+l 

ij 


+ d 


n.n+l 

i-lj 


+ e 


n,  n+1 


ij-i 


(3-9) 


where  an,  bn,  c11  etc.  are  coefficient  matrices  containing  only  n-level 
quantities,  fn  is  the  finite  difference  form  of  the  explicit  part  of 
equation  (3-1),  and  <j>^  is  the  column  vector  containing  the  dependent 
variables  at  point  ij . If  a single  row  or  column  of  grid  points  were 
being  solved  by  equation  (3-9)  (as  in  a one-dimensional  problem)  the 
result  would  be  a block-tridiagonal  matrix  which  could  be  quickly  solved 
by  a standard  elimination  technique.  Application  of  equation  (3-9)  to 
a field  of  many  rows  or  columns  of  points  results  in  a large  cumbersome 
matrix  which  can  be  solved  by  Gauss  elimination  or  some  iterative  tech- 
nique. The  computation  time  required  for  solution  by  either  method 
increases  rapidly  with  the  size  of  the  grid. 

As  mentioned  previously,  the  computation  time  of  solution  for  a 

one-dimensional  problem  is  retained  by  use  of  ADI  techniques . Of  the 

31 

many  ADI  schemes  (see,  for  example,  Yanenko  ),  Briley  and  McDonald  use 
a form  of  the  general  procedure  of  Douglas  and  Gunn.  This  procedure 


l 
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generates  the  ADI  scheme  as  perturbations  of  the  fundamental  implicit 
scheme.  It  is  a multistep  method  (one  step  for  each  spatial  dimension) 
where  the  first  step  approximates  the  implicit  equation  and  subsequent 
steps  add  corrections.  Yanenko  calls  this  the  "method  of  stabilizing 
corrections"  and  shows  that  the  method  has  the  two  important  properties 
of  consistency  and  stability.  Briley  and  McDonald  point  out  that  the 
consistency  property  allows  for  the  use  of  physical  boundary  conditions 
for  the  intermediate  step  with  no  loss  in  accuracy  for  steady  state 
solutions. 

Each  step  of  the  procedure  involves  the  implicit  solution  in  one 
of  the  coordinate  directions.  This  results  in  a system  of  one-dimen- 
sional, block  tridiagonal  matrices  which  are  easily  solved  by  standard 
block  elimination  methods.  As  an  example,  applying  the  ADI  scheme  to 
equation  (3-1)  gives: 


_ , * n , n * nn,  ..  r n n, 

-6x[p  u+pu-pu]^-  5y[p  v ] ^ 


(3-10) 


** 

(p “ 


. , * n , n * n n,  . r **  n . n **  n n-, 

<5x[p  u +p  u -p  u J^j-SyiP  v +p  v -p  v Jjj  (3-11) 


where  * Indicates  the  intermediate  quantity  given  by  equation  (3-10)  and 
**  represents  the  quantity  evaluated  by  equation  (3-11).  The  complete 
set  of  ADI  equations  is  given  in  Appendix  C.  Now,  equations  (3-10)  and 
(3-11)  can  be  written  as 


Dxn(j)*  + Dyn<f>n  + Sn 


(3-12) 
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**  n 

„0  — o v n * n ** 

(a — rr-M  - D*  ♦ + Dy  <j>  + Sn  (3-13) 

AC  ±j 

where 

Dx11  - |-  6x{un}  -<Sx{pn}  0 | 

Dy"  - |-  5y{vn}  0 -6y{pn}| 


Sn  ■ <Sx(pnun)  + 5y(pnvn)  <j>  = |^| 

and 

{ } implies  multiplication  with  $ before  the  difference  opera- 
tion is  applied. 

£ ick 

Now  <j>  is  computed  in  the  intermediate  equation  (3-12)  and  $ 

27  ** 

is  obtained  in  (3-13).  According  to  Douglas  and  Gunn.  ./p  is  within 
Q(At^)  of  4>n  ^ and  so  is  taken  at  <J>n+^.  The  system  can  be  simplified 
by  subtracting  (3-12)  from  (3-13)  to  get  the  new  system: 


<a-irE->u  ' D*n»ij + 

+ sn 

(3-14) 

**  * 

- *«> 

(3-15) 

Equations  (3-14)  and  (3-15)  can  be  written  as: 

— n * =n  * — n * 

Ci  ^i-lj  + bi  *ij  + ai  ^i+lj 

(3-16) 

= n **  = n **  =n  ** 

— * — n 

(3-17) 

Yj  ♦ ij-l  + Si  *ij  + ai  *ij+l  “ 

where  c°  , b ",  etc.  are  the  coefficient  matrices  of  the  unknowns 
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— n — * 


— n 


$>„,  etc.  The  column  vectors  d^  , , and  ru  contain  the 

explicit  terms.  The  order  of  the  coefficient  matrices  is  equal  to  the 
number  of  dependent  variables.  The  forms  of  a”  , b^1  , c^1  , etc.  come 
from  the  finite  difference  equations,  and,  for  the  interior  points,  all 
have  the  same  form.  (See  Appendix  D.)  For  points  adjacent  to  a 
boundary,  the  finite  differencing  must  be  modified.  These  modifications 
are  discussed  in  Chapter  V. 

The  solution  for  a single  time  step,  then,  proceeds  as  follows: 

1.  Equation  (3-16)  is  applied  at  successive  rows  (x-direction) 

to  generate  a series  of  coupled,  one-dimensional  equations 

(there  being  one  set  of  coefficient  matrices  a”  , , c^1 

and  a d”  for  each  point  in  the  row) , which  are  arranged  into 

one  block-tridiagonal  matrix  for  each  row.  The  matrix  is 

then  solved  by  a standard  block  elimination  technique  (see 

* 


Chapter  IV)  to  give  the  values  of  . 

2.  The  second  step  is  similar  to  the  first  except  that  equation 
(3-17)  is  applied  to  successive  columns  (y-direction) , which 


gives  the  <j>  vector  for  the  flow  field. 

It  should  be  noted  that  the  "splitting"  of  the  Douglas-Gunn  pro- 
cedure can  be  done  in  any  coordinate  direction  and  does  not  require 
association  with  coordinate  directions.  The  criteria  used  is  that  the 
associated  matrices  are  easily  solved.  The  mixed  derivatives  can  be 
treated  implicitly,  therefore,  but  this  Increases  the  number  of  inter- 
mediate steps  and  greatly  complicates  the  procedure  and  was  not  done 
here.  Thus,  as  the  method  was  applied,  it  was  not  totally  implicit. 
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Experience  has  shown  that  this  does  not  seem  to  hinder  the  ability  of 

the  scheme  to  use  larger  time  steps  than  those  required  by  the  viscous 

stability  requirement.  Briley  and  McDonald,  for  example  were  able  to 

2 

use  At  = 20.6  Aty  (where  Aty  = (Ax)  /2v)  and  At  = 1471  At  . They 

Cr  L 

point  out  that  explicit  treatment  of  a dissipation  term  and  V-u  appeared 
not  to  affect  the  stability  of  the  procedure,  even  at  these  large  time 
steps.  Hence,  it  was  deemed  not  necessary  to  compute  mixed  derivatives 
implicitly.  Subsequent  experience  with  computations  for  supersonic 
base  flow  confirmed  this. 
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CHAPTER  IV 

MATRIX  SOLUTION 

Equations  (3-16)  and  (3-17)  generate  a series  of  coupled,  linear, 

one-dimensional  equations  for  each  row  and  each  column  of  grid  points  in 

the  flow  field.  Each  series  of  equations  represents  a complete  block- 

tridiagonal  matrix  system  which,  as  previously  mentioned,  can  be  solved 

by  standard  techniques.  The  one  used  here  is  the  L-U  decomposition  and 

28 

back-substitution  (LUBS)  method  described  by  Isaacson  and  Keller. 

/ 

To  illustrate  the  procedure  consider  the  solution  along  the  jth 
row  of  grid  points.  At  each  point,  Eq.  (3-16)  gives: 

=n  * , =n  * , =n  * -rn 

Ci  *i-lj  + bi  *ij  + ai  *i+lj  = di 

For  the  whole  row,  the  system  emerges  as: 


, n n 

,n  i 

bi  *i 

♦ij 

di 

n , n n 

* 

,n 

C2  b2  a2 

*2j 

d2 

n . n n 

* 

,n 

C3  b3  a3 

• 

• • 

*3J 

• 

d3 

• 

■ • 

. * n 

• 

• 

* 

,n 

“N-l 

S-l 

n , n 

* 

,n 

% bN 

* Nj 

Sr 
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where  N is  the  number  of  grid  points  in  a row.  Now  each  of  the  coeffi- 
n , n n 

cients  Ci  ’ » aad  are  themselves  square  matrices  of  the  order 

k,  the  number  of  dependent  variables. 

The  first  step  is  to  convert  the  coefficient  matrix  to  the 
product  of  the  upper  and  lower  triangular  matrices  (hence  the  name  LU 
decomposition) . 


where  I “ identity  matrix. 

This  is  accomplished  by  using  the  recursion  formulas: 


bi-  bl 

al = bl  al 

(4-1) 

bi  " bi  * ciaI-i 

(i  - 2,3 N) 

(4-2) 

ai  - <b;>-\ 

(i  - 2,3, .. . ,N-1) 

(4-3) 

Hence,  the  system  becomes: 
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k 


- 

L 


which  is  solved  as 


Uy  - d ; L<{>  = y ; y - 


(4-5) 


Recursions,  for  y and  $ are: 


yx  - <b^)~1d1;  y±  - (b’)"-L(di  - c.y^)  (i  - 2,3,. ..,N)  (4-6) 


v-1 


and 


“ yN’  *1  " yi  " ai$i+i  (i  * N-l,  N-2,...,l)  (4-7) 

This  method  is  seen  to  take  full  advantage  of  the  large  number  of  zeros 
in  the  matrix  by  performing  operations  only  on  the  nonzero  elements  of 
the  coefficient  matrix.  It  is  thus  seen  as  particularly  efficient  and 
suitable  for  use  in  a computer  program.  Further  computational  efficiency 
is  gained  if  equations  (4-1),  (4-3)  and  (4-6)  are  solved,  not  as  written 
by  inverting  the  b^  matrix,  but  by  solving  the  linear  system  with  b^,  on 
the  left  hand  side  of  the  equation.  Gauss  elimination,  for  example, was 
used  in  this  problem,  though  other  techniques  could  be  used.  The  ques- 
tion arises,  in  the  solution  of  the  1 linear  system  with  on  the 
left  hand  side,  as  to  whether  a pivoting  strategy  would  aid  in  the 
reduction  of  any  round  off  error.  These  round  off  errors  may  arise,  for 
example,  from  the  fact  that  so  many  arithmetic  operations  are  being 
performed  even  in  the  solution  of  one  row  of  grid  points.  A check  for 
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round-off  errors  was  made  for  this  problem  by  making  a single  computa- 
tion using  double  precision  (single  precision  on  the  CDC  6600  of  the 
Georgia  Tech  CYBER  74  is  14  decimal  places).  This  computation  was  run 
for  100  time  steps  and  the  results  were  identical  to  a computation  made 
using  single  precision  arithmetic.  Thus  no  pivoting  appeared  necessary. 
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CHAPTER  V 

BOUNDARY  AND  INITIAL  CONDITIONS 

The  finite  difference  equations  all  have  the  same  form  when 
applied  to  interior  grid  points.  When  a cell  is  adjacent  to  a boundary 
these  equations  must  be  modified.  This  chapter  presents  the  modified 
forms  of  the  finite  difference  equations  for  the  various  boundaries 
encountered  in  the  flow  field. 

While  all  flows  are  governed  by  the  same  set  of  equations,  the 
variety  of  phenomena  (bubbles,  shocks,  recirculation,  etc.)  arise  due 
to  differing  boundary  conditions.  It  follows  that  these  boundary  con- 
ditions in  a numerical  study  must  be  specified  carefully,  as  was  indeed 
discovered.  The  cell  integration  formulation  affords  great  conceptual 
aid  here,  in  that  it  becomes  clear  which  group  of  terms  needs  to  be 
modified.  The  problem  becomes  one  of  choosing  among  several  plausible 
forms.  Most  of  these  forms  are  outlined  in  Roache. 

Complicating  the  matter  is  the  fact  that  appropriate  forms  appear 
to  vary  with  solution  procedure.  Roache  cites  several  instances  where 
one  form  gave  good  results  for  some  methods,  but  caused  numerical  diver- 
gence in  others  (see,  for  example,  Roache, p.  280).  Allen  and  Cheng 
even  found  that  near-wall  flux  terms  needed  to  be  modified  within  the 
same  method  for  a finer  mesh  to  get  physically  meaningful  densities. 

This  is  not  to  imply  that  only  a single  form  will  work  in  a given  condi- 
tion; just  that  the  solution  can  be  very  sensitive  to  form  change.  The 
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solution  can  also  be  quite  insensitive  to  different  forms.  Many  forms 
for  the  viscous  terms  near  and  on  the  walls  were  tired,  and  gave  negli- 
gible changes  in  the  solution.  It  is  not  surprising  then,  that  a major 
effort  was  required  to  arrive  at  appropriate  modifications  of  the  finite 
difference  equations  for  the  edge  regions. 

The  remainder  of  this  chapter  presents  the  final  forms  of  the 
boundary  conditions  used  for  solution.  A more  complete  discussion  of 
other  boundary  condition  forms  considered  is  given  in  Chapter  VI, 
Results.  Only  the  conditions  for  laminar  flow  are  given  here.  The 
extension  of  this  method  to  turbulent  flow  and  associated  boundary  con- 
ditions is  given  in  Bangert  and  Roach. 


5.1  Upper  Wall 

The  upper  wall,  labeled  BC  in  Figure  3,  is  no-slip,  impermeable, 
and  adiabatic.  A cell  adjacent  to  BC  (see  Figure  5)  must  then  have 


(u)  - 0 

y- 


(v)y.  - 0 


Oe/3y)y_  = 0 


(5-1) 


(5-2) 


(5-3) 


Equations  (5-1)  and  (5-2)  imply: 


(3u/3x)y_  - 0 


(3v/3x)y_  - 0 


(5-4) 


(5-5) 


Upper  wall 
along  BC 


i,j+l 

• 

• 

y+ 

X-  xl 

y- 

1+1.  J 
• 

- 

i.j-i 

i+1.3-1 

• 

Back  wall 
along  CD 


Figure  5.  Boundary  Cells  - Upper  Wall,  Back  Wall  and  CenCerline 
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Similarly,  the  shear  work  terms  (from  the  energy  equation)  on  the  wall 


are  zero: 


[■£-  (3u/3y  + 3v/3x)u]  = 0 

Re  y— 


(5-6) 


(y  3v/3y  - y 3u/3x)v]y_  * 0 


(5-7) 


The  nondimensional  pressure  at  the  wall,  (pe)^_ , was  evaluated 
by  a linear  extrapolation  through  (ij)  and  (i,j+l)  giving. 


(pe)y-  * 2[3(pe)ij  ~ (pe)ij+l] 


(5-8) 


This  technique  was  also  used  by  Allen  and  Cheng. 

When  normal  derivatives  are  required  at  the  surface,  a second- 
order  accurate,  one-sided  difference  was  used: 


<3*/3>V  - 3^7  (-8V  + 9*ij  “W 


(5-9) 


(— ) 

Wy- 


3Ay  (9uij  " Uij+1) 


(5-10) 


(— ) 

Wv~ 


3Ay(9vij  ' Vij+1) 


(5-11) 


5.2  Back  Wall 

This  wall,  labeled  CD  in  Figure  3,  is  also  Impermeable,  no-slip, 
and  adiabatic.  A cell,  then,  with  the  wall  on  the  x-edge  (see  Figure  5) 
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(u)  = 0 

x- 

(v)  = 0 

X- 

and, 

Oe/3x)x_  = 0 


(5-12) 

(5-13) 

(5-14) 


Equations  (5-12)  and  (5-13)  imply  that: 


(3u/3y)x_  =*  0 (5-15) 

Ov/3y)x_  = 0 (5-16) 


also: 


r V A 3u 
lRe  v3  3x 


! I?)u]x-  - 0 


(5-17) 


rJL.  (iH  + lZ)v] 
lRe  3y  3x'  Jx- 


(5-18) 


using  the  second-order  form  for  the  normal  derivatives,  similar  to 
Equation  (5-9),  gives: 


and 


(3u/3x)x 


Ui+ij ) 


(5-19) 


<3v/3x)x-  ' <9vlJ  ~ vi+ij) 


(5-20) 


The  form  used  for  the  pressure  on  the  back  wall  was  not  the 
extrapolation  used  as  on  BC,  but  rather: 

(pe)x_  - (9pe^j  - pe±+lj) 


(5-21) 
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Note  that  this  form  is  equivalent  to  a zero  pressure  gradient  normal  to 
the  wall  at  x-using  the  second-order  form  for  normal  derivatives.  The 
reason  for  the  use  of  this  form  is  given  in  Chapter  VI  in  the  section 
on  boundary  conditions. 


5.3  Centerline 

The  centerline,  labeled  DE  in  Figure  3,  is  a plane  of  symmetry 
and  thus  has  no  mass  flux  across  it.  A typical  cell  (see  Figure  5)  is 
adjacent  to  DE  at  its  y-  edge.  So 


(v)y_  - 0 (5-22) 

(3v/3x)y_  = 0 (5-23) 

and 

(3<J>/3y)  - 0 (4>  ^ v)  (5-24) 


The  second  order  form  for  the  derivatives,  eq.  (5-9)  is  used 
for  the  normal  derivative  of  v: 


Ov/3y)y_  * 3^  - vlj+1> 


(5-25) 


and  also  gives  a consistent  form  for  the  nonzero  variables  at  y-: 

N <*V  - 1 <5-26> 

Equations  (5-22)  through  (5-24)  imply  that  the  shear  work  terms 
adjacent  to  DE  are  also  zero: 

["Re  (3u/3y  + 3v/3x)u]y_  * 0 


(5-27) 
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(-|  3v/3y  - ~ 3u/3x)v]y_  = 0 


(5-28) 


5.4  Outflow  Boundary 

The  outflow  boundary,  EF  in  Figure  3,  is  adjacent  to  the  outflow 

cells  at  their  x+  edges.  Little  can  be  assumed  about  the  flow  here 

regarding  specification  of  gradients  or  variables  because  the  flow  is 
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not  known  a priori.  For  explicit  schemes,  Roache  suggests  various 
extrapolation  methods.  He  notes  that,  in  most  cases,  linear  extrapola- 
tion is  satisfactory,  except  perhaps  when  a shock  crosses  the  boundary. 
Thus,  an  implicit,  linear  extrapolation  procedure  was  used  to  specify 
conditions  at  the  x+  edge  of  the  outflow  cell.  The  schemes  which  com- 
puted the  conditions  at  the  outflow  cell  points  by  explicit  extrapola- 
tion all  caused  divergence  for  At  > At„„  . Zero  gradient  forms  for  the 
outflow  cell  points,  both  explicit  and  implicit,  resulted  in  wiggles  in 
the  steady  state  solution.  Details  are  described  in  the  section  on  bound- 
ary conditions  in  Chapter  VI. 

The  values  of  all  the  dependent  variables  at  x+  are  obtained  by 
a linear  extrapolation  from  ij  and  i-l,j: 


*x+  “ 2 (3*ij  ' *i-l,J‘ 


(5-29) 


Thus  normal  derivatives  at  x+  are  equal  to  the  normal  derivative  at  x-: 


W/5,)*.  - <»♦/»*>*.  - H <♦«  - (5-30> 


5.5  Upper  Boundary 

For  the  upper  boundary,  labeled  AF  in  Figure  3,  the  properties 
were  determined  explicitly  (i.e.,  after  the  rest  of  the  flow  field  has 
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been  computed)  by  the  simple  wave  procedure  used  by  Allen  and  Cheng 

and  outlined  by  Roache  (Reference  30,  pp.  282-283).  Briefly,  it  was 

assumed  that  properties  are  constant  along  the  straight,  left-running 

characteristic  line  passing  through  an  upper  boundary  point  ij . This 

line  is  determined  by  the  angle  (p  +9),  where  p * arc  sin (1/M)  is 

in  m 

the  local  Mach  angle  (M  ^ 1)  and  9 - arc  tan(v/u)  is  the  local  flow 
direction.  The  properties  on  the  characteristic  line  are  determined  by 
linear  interpolation  between  (i-1,  j-1)  and  (i,j-l)  or  (i-l,j-l)  and 
(i-l,j)  depending  on  the  local  (pm  + 9)  and  the  ratio  Ay/Ax.  Figure  6, 
shows  the  two  cases. 

For  tan  (p  + 9).  . . . > Ay/Ax,  the  characteristic  line  runs 

m 1— 1, j— 1 

between  (i-1, j-1)  and  (i,j-l),  and  the  properties  at  P are  determined 
by: 

*P  " * i-i, j-i + (-ir)(*i,j-i " +i-i,j-i)  (5_31) 


Then  * <jip  is  used  to  solve  for  the  upper  boundary  points.  The  value 
of  l and  thus  the  location  of  P are  determined  as  follows.  Consider  the 
quantity 


u - tan  [90°  - (p  +9)] 
m 

_ _ Ax- l 

By  geometry:  Up  “ 

But  from  (5-40) 


<D 

P 


“i-1, j-1  + (Ax)(“i,j-l  ~ “i-lj-l5 


(5-32) 


(5-33) 


Figure  6.  Boundary  Cells  - Upper  Boundary. 
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Equating  (5-41)  to  (5-42)  and  solving  for  gives. 


(Ax/Ay  - 

to  K.j-l  ' “i-l.J-15  + ^ 


(5-34) 


For  the  case  of  tan  (ym  + 9)  < ^ , the  characteristic  line  runs 


m 


between  (i-l,j-l)  and  The  properties  at  p are  determined  by 


’i-l.j-l  + Ay  (*i-l,j  " ♦i-l.J-l5 


(5-35) 


and  then  <t>  . =*  $ as  before.  Here  the  quantity  w = tan  (y  + 9)  is  used 
lj  p m 

and 


is  equated  with 


u 


P 


Ay  - l 
Ax 


“i-i.j-1  + a7  (“i-ij 


“i-ltJ-l) 


to  get 


(Ay/Ax  - 

Ay  (“i-l,j  “ “i-l.j-l5  + Ax 


(5-36) 


(5-37) 


(5-38) 


5.6  Inflow  Boundary 

The  flow  properties  on  the  inflow  boundary  (AB  in  Figure  3) 
were  held  fixed  during  the  computation.  Their  values  were  chosen 
starting  with  the  assumed  velocity  profile  used  by  Allen  and  Cheng. 


u « j (2n7  - 7n^  + 14o)  0 <_  n <_  1 

u - 1 n > 1 


(5-39) 
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where  n ■ y/6 

6 ■ nondimens ional  boundary  layer  height 
The  internal  energy  profile  was  determined  from  the  Busemaun  integral 
of  the  compressible  laminar  boundary-layer  energy  equation  for  an  adia- 
batic wall: 


e - 1 + 7 (y  - DM,2  (1  - u2) 


(5-40) 


with  the  assumption  of  constant  pressure  through  the  inflcw  boundary 
layer,  the  density  is  determined  from  the  equation  of  state.  Here 


p “ 7 


(5-41) 


With  these  values  fixed  the  vertical  velocity  component  is  determined 
from  the  boundary  layer  equations.  As  in  Allen  and  Cheng,  an  ordinary 
differential  equation  for  v can  be  derived  by  combining  the  x-momentum 
and  continuity  equations  and  using  the  energy  integral.  Integrating 
both  sides  gives: 


2 

v(n)  * - \ [1  +y(y-l)M12  (1-hi2)]  i-|  dn  0<n<l 

« u dn 


(5-42) 


v(n)  - v(l) 


where  q ■ y/6  5 ■ .41 


n > 1 


5.7  Initial  Conditions 


The  unsteady  equations  require  that  initial  conditions  be  speci- 


fied everywhere  before  computation  can  begin.  This  specification  is 


arbitrary  but  some  care  must  be  taken.  Initial  conditions  with  very 
steep  gradients  near  the  expansion  corner,  for  example,  were  found  to 
cause  divergence.  An  examination  of  the  effect  of  a few  different 
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initial  conditions  is  discussed  in  Section  6.3-4.  For  the  development 
of  the  procedure  and  for  most  of  the  computations  the  same  initial  con- 
ditions were  used.  The  boundary  layer  and  freestream  conditions  were 
imposed  along  the  upper  wall  to  corner  point  C.  Beyond  this  corner  and 
above  it,  the  freestream  conditions  were  applied  but  with  v ■ 0.  Below 
the  expansion  corner  v was  zero  also  and  u was  30Z  of  the  freestream 
value.  The  parameters  of  the  flew  were  set  to  correspond  to  Allen  and 
Cheng's  Case  Bl,  where  the  freestream  Mach  number  was  3,  the  Reynolds 
number  was  550,  the  nondimensional  boundary  layer  height  was  0.41,  and 
the  walls  were  adiabatic.  This  also  corresponds  to  Allen's  cases  Cl 
and  C5. 


L 
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CHAPTER  VI 

RESULTS 

This  chapter  is  divided  into  three  parts.  The  first  part  discusses 
the  various  finite  difference  schemes  tried  for  the  boundary  conditions 
and  the  reasons  for  the  final  choices.  The  second  part  gives  the 
results  of  the  flow  field  computations  and  the  comparison  with  the 
results  obtained  by  Allen  and  Cheng.  The  last  part  discusses  the 
results  of  some  numerical  tests  on  the  method,  especially  the  use  of  a 
finer  grid,  larger  time  steps,  and  different  initial  conditions. 

6.1  Boundary  Conditions 

Like  the  finite  difference  schemes,  the  boundary  conditions  can 
be  either  explicit  or  implicit.  For  an  overall  implicit  procedure,  it 
is  desirable  to  have  implicit  boundary  conditions  to  prevent  problems 
associated  with  time  lagging  of  the  boundary  conditions  behind  the  flow 
field.  For  time  steps  much  larger  than  the  explicit  stability  limit, 
the  use  of  explicit  forms  on  some  boundaries  may  not  be  possible.  Addi- 
tionally, the  use  of  implicit  boundary  conditions  may  accelerate  conver- 
gence, a very  desirable  feature.  When  the  term  "implicit"  is  applied  to 
any  of  the  forms  to  be  described  below,  it  means  that  the  form  was 
incorporated  into  the  block  tridiagonal  matrix  and  the  variables  at  the 
boundary  solved  along  with  those  on  the  interior.  "Explicit"  refers  to 
those  schemes  which  compute  the  boundary  variables  after  those  at  inter- 
ior points  have  been  computed. 
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As  stated  in  the  previous  chapter,  boundary  conditions  can  not 
be  treated  haphazardly.  Care  must  be  taken  in  the  selection,  from  the 
many  forms  possible,  to  prevent  "wiggles"  or  even  divergence.  Several 


of  these  forms  are  discussed  in  the  following  sections  but  some  prelimi- 
nary comments  about  wiggles  are  appropriate. 

Wiggles  are  nonphysical  spatial  oscillations  occurring  in  the 
solution.  Roache  points  out  that  wiggles  are  not  usually  caused  by 
iterative  instability,  nonlinearities,  or  spatially  varying  coefficients, 
but  are  actually  the  solution  of  the  finite  difference  equations. 

Moretti  gives  several  examples  where  the  appearance  of  wiggles  was 
caused  by  poor  modeling  of  the  physical  behavior  (or,  in  some  cases,  no 
modeling  at  all),  particularly  poor  treatment  of  the  boundary  conditions. 

Standard  techniques  for  treating  procedures  which  produce  wiggles 
are  the  use  of  artificial  viscosity  or  switching  from  central  differ- 
encing to  upwind  differencing.  Roache,  however,  shows  that  the  two  are 
roughly  equivalent,  the  truncation  error  of  upwind  differencing  corre- 
sponding to  an  artificial  viscosity.  Use  of  upwind  differencing  can 
become  complicated  in  regions  where  the  flow  direction  is  not  known 
beforehand,  as  in  a recirculation  region.  Further,  artificial  viscosity 

may  tend  to  smooth  important  flow  features  over  too  large  an  area  or 
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even  change  the  flow  problem.  Moretti,  for  example,  shows  that  artifi- 
cial viscosity  completely  wipes  out  a shock  in  a Laval  nozzle  and  causes 
the  procedure  to  be  significantly  in  error  for  the  critical  throat  con- 
ditions. 


Because  of  the  foregoing,  the  assumption  was  made  that  appropriate 
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forms  for  the  boundary  conditions  could  be  found  without  having  to 
abandon  central  differencing  or  having  to  apply  artificial  viscosity. 

This  proved  to  be  the  case.  Indeed,  all  nonphysical  behavior  of  pre- 
liminary solutions  came  from  ill  treatment  of  one  or  more  of  the  bound- 
ary conditions.  (Many  times  a programming  error  connected  with  the 
boundaries  was  responsible.  A considerable  amount  of  time  was  required 
to  completely  "debug"  the  program,  due  largely  to  the  complexity  of  the 
equations) . 

As  mentioned  before,  a great  advantage  in  the  cell  integration 
technique  is  the  conceptual  aid  afforded  in  specifying  boundary  condi- 
tions. When  fluxes  are  required  to  be  zero,  for  example,  it  is  clear 
which  terms  should  be  eliminated  from  the  finite  difference  equations 
so  that  enforcement  is  automatic.  Ambiguity  arises  when  a dependent 
variable  or  its  derivatives  are  not  specified  on  the  boundary  by  physi- 
cal conditions  (such  as  wall  pressure).  Then,  one  of  the  many  extrapo- 
lation or  one-sided  forms  must  be  chosen.  As  stated  in  Chapter  V,  much 
effort  was  required  to  determine  which  forms  were  stable,  consistent, 
and  accurate.  This  section  summarizes  that  effort . 

6.1-1  Wall  Boundaries 

Surprisingly,  the  solution  was  insensitive  to  many  of  the  differ- 
ing forms  for  conditions  at  the  walls.  Whenever  a problem  occurred,  it 
was  either  a programming  error  or  a problem  on  some  other  boundary.  This 
was  perhaps  a result  of  no  slip  and  impermeability  being  strictly 
enfroced.  These  are  the  dominant  conditions  that  characterize  the  wall. 

Most  of  the  forms  used  on  the  wall  were  those  used  by  Allen  and 
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Cheng,^  although  many  other  forms  were  tried  attempting  to  get  rid  of 
wiggles  in  the  solution  near  the  back  wall.  It  was  discovered,  however, 
that  these  wiggles  were  a result  of  poor  downstream  boundary  treatment 
(see  Section  6.1-3)  and  the  Allen  and  Cheng  forms  appeared  to  work  best. 
The  different  forms  used  seemed  to  have  little  effect  on  the  solution 
although  second-order  forms  gave  more  physically  correct  answers  than 
first  order  forms,  as  is  explained  next. 

Since  the  walls  are  impermeable,  there  is  no  mass,  momentum,  or 
energy  flux  across  a cell  boundary  coinciding  with  a wall.  Hence,  the 
contribution  to  the  convection  term  from  the  cell  edge  is  zero,  leaving 
a first  order,  one-sided  difference  for  the  flux  derivative.  For  exam- 
ple, consider  the  finite  difference  form  for  a flux  tern  of  a cell  ij 

mass 

x-mom . 

(6-1) 

y-mom. 
energy 

For  a wall  at  x-,  (i.e.,  ij  is  next  to  the  back  wall)  (pu$)x_  - 0 and 

(3pu*/3x)  = i (pu^  - 0)  > Kpub)iuJ  + (6-2) 

When  this  form  was  used  the  zero  streamline  extended  from  the 
rear  stagnation  point  forward  to  another  point  on  the  centerline.  The 
expected  behavior,  and  that  obtained  by  Allen,  Allen  and  Cheng,  and 
Kronzon,  et  al.,  was  for  the  zero  streamline  to  extend  forward  from  the 
rear  stagnation  point  to  a point  on  the  back  wall  just  below  the  expansion 


(3pu$/3x)^  = (pu<J>v+  - pu$x_)  6 
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corner.  The  reason  for  the  unexpected  streamline  behavior  can  be  seen 
in  the  velocity  vector  plot  in  Figure  7.  As  will  be  explained  again  in 
Section  6.2,  the  plot  shows  total  velocity  magnitudes  and  directions 
(with  the  arrow  heads  a constant  size)  for  each  point  in  the  coarse  mesh 
field.  The  zero  streamline  intersects  the  centerline  near  point  D 
instead  of  on  the  wall  near  C because  all  the  velocities  in  the  column 
nearest  the  back  wall  are  positive,  when  most,  starting  at  the  point  near- 
est D should  have  been  negative.  That  equation  (6-2)  was  the  problem  can 
be  seen  by  considering  a typical  flux  term  for  a cell  adjacent  to  the 
back  wall. 


h (pu*}ij  "ir  I(pu^x+ 

= 't  (pU*}x+ 


(?u?)x_] 


(6-3) 


since 


(pud)x_  * (PUl^waii  ” 0 

Using  equation  (6-2)  for  (pu$)x+  assumes  a linear  variation  in  the 
flux  across  the  cell.  Thus 

PU(()  - CjX  - j (P^i+ij  + ou?ij)x  (6-4> 

where  x - distance  from  the  wall. 

This  means  that 


3x  (pU*> 


I (9U*i+lj  + pu*iJ> 


(6-5) 
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Figure  7.  Velocity  Vectors  for  First-Order  Backvall  Flux  Differencing. 


45 


for  the  whole  cell  from  the  wall  to  x+.  From  the  steady-state  contin- 
uity equation,  evaluated  at  the  wall 


(^~)  + (*£)  - 0 


3x 


w 


but 


■ o 


since  pv  - 0 all  along  the  wall.  So 


(|^)  - 0 

3X  w 


Now 


. u l±+.i£”i 

3x  p 3x  v3x 


Evaluating  equation  (6-8)  at  the  wall  gives 


<i£M,  .0 

3X  W 


since  (u)  ■ 0 and  0§^)  = 0.  But  this  implies  that 


3x 


w 


and 


C1  * J + pu<,>ij^  * 0 


pu*i+lj  " * pu*ij 


(6-6) 


(6-7) 


(6-8) 


(6-9) 


(6-10) 


(6-11) 


Since  velocities  are  small  near  the  lower  part  of  the  back  wall,  the 
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dominant  flux  terms  are  those  for  <j>  = 1 and  <j>  ■ e,  the  mass  and  energy 
fluxes.  Therefore,  since  p,  1,  and  e are  all  positive,  u is  required 
to  change  sign  between  (i,j)  and  (i+l,j).  Thus  u was  positive  for  the 
column  nearest  the  wall  and  negative  in  the  next  for  the  cells  in  the 
recirculation  region.  Negative  u was  expected  in  both.  Further,  the 
fact  that  pu<|>  changes  sign  between  (i,j)  and  (i+1, j ) , having  been  zero 
at  the  x-edge  of  (i,j),  implies  at  least  quadratic  behavior  in  violation 
of  the  linear  assumption  of  equation  (6-2). 

Allen  and  Cheng  encountered  negative  densities  on  the  upper  part 
of  the  back  wall  near  the  expansion  corner  when  they  used  equation  (6-2) 
and  fine  grid  spacing.  Equation  (6-11)  may  explain  this  behavior.  In 
this  region,  u is  not  small  so  that  for  =»  u,  equation  (6-11)  implies 

(pu2)±j  - -(pu2)±+lj 

which  may  explain  why  negative  densities  were  encountered.  It  has 
already  been  shown  that  the  assumption  of  a linear  variation  in  pu<j>  leads 
to  a violation  of  that  assumption  and  some  physically  unrealistic  results. 
Equation  (6-9)  and  the  fact  that  pu<ji  ■ 0 at  the  wall  imply  at  least  a 
quadratic  variation  in  the  flux  across  the  cell.  Allen  and  Cheng  sug- 
gested a second-degree  polynomial  through  the  two  cell  centers  nearest 
the  wall,  and  used  the  known  value  of  u ■ 0 at  the  wall  to  get  an  expres- 
sion for  the  flux  at  x+: 

(pu*>x+  " 7 [(pu^iarfl,j  + 3(pU<,,)ia.,j]  (6'12) 


This  gives 
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Opu*/3x)i(i)>j 


3Ax  ««”*> 


iu+1,  j 


+ 3<pu*)i»,jI 


(6-13) 


This  form  gave  the  correct  behavior  for  the  cells  adjacent  to  the  back 
wall  in  the  present  calculations. 

The  pressure  at  the  wall  is  unknown.  Allen  and  Cheng  suggested 
a linear  extrapolation  from  the  two  nearest  cells: 


(6-14) 


This  form  was  used  initially  on  the  back  wall  and  its  equivalent  form 
was  used  on  the  upper  wall. 

Two  other  forms  were  tried.  Allen  and  Cheng's  results  showed 
that  the  pressure  gradient  normal  to  the  back  wall  is  very  small.  Thus, 
the  other  two  forms  tried  for  the  wall  pressure  involved  a zero  pressure 
gradient.  The  first  form  is  the  simple 


The  second  form  makes  use  of  the  second-order,  one-sided  form  (equation 
(5-21)) 


<0e>u, 


1 

8 


[9(pe) 


iw,j 


(pe)iw+l,j] 


(6-16) 


Experience  with  calculations  on  all  the  boundaries  indicates  that  a sec- 
ond-order form  is  superior  to  a first-order  form  in  that  the  results  are 
usually  more  physically  realistic  and  more  consistent  with  the  basic 
finite  difference  method.  Indeed,  when  equation  (6-2)  was  used  for  the 

i 

flux  terms,  equation  (6-16)  tended  to  bring  the  zero  streamline  off  the 
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centerline  and  back  to  the  back  wall.  Equations  (6-14)  and  (6-15),  on 
the  other  hand,  resulted  in  little  change.  When  equation  (6-12)  was 
used  for  the  flux  terms,  there  was  little  difference  in  the  solution  with 
any  of  the  wall  pressure  forms.  Thus,  for  some  conditions  equation 
(6-16)  gave  better  physical  results  and  was  the  form  used  for  the  back 
wall  pressure. 

Four  forms  were  tried  for  the  viscous  terms  involving  cross-deriva- 
tives. None  of  these  forms  appeared  to  be  preferred.  Most,  in  fact, 
gave  identical  answers  to  several  decimal  places.  There  are  two  types 
of  cross-derivatives.  Considering  the  upper  wall,  the  first  form  is 


3y 


(£  where  $ 

3x 


£ 


u/Re  momentum  eqns. 
(u/Re)$  energy  eqn. 


(6-17) 


Applying  cell  integration  this  becomes 


_1_ 

Ay 


-=-{  (-  12.)  ] 


7+ 


(6-18) 


since 


(li) 

y_ 


0 


Two  forms  for  equation  (6-18)  are: 


_i_  (r  M.)  _ -i-(c)  (2£.)  . ^ 1 f1- * ? J LJ  — iniJ)] 

Ay  U d*  y+  Ayvuy+Sx'y+  Ay  ^ 2 n 4ax 


(6-19) 
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-L(rii)  + r /li+illiizii)!  (6_20) 

AyU  3x'y+  Ay  2Kij+lV  2Ax  ; Hj 1 2Ax  ( } 

The  question  here  is  whether  to  split  (£  -j^-)  , into  (5)  , Cr^O  , or  to 

3x  y+  y+  3x  y+ 

evaluate  it  as  a single  quantity.  Naturally,  the  use  of  either  (6-19) 
or  (6-20)  depends  on  the  form  used  for  the  interior  points,  to  be  con- 
sistent. Both  forms  were  tried  and  the  results  were  identical.  Equa- 
tion (6-20)  was  the  final  form  used  because  the  non-cross-derivative 
viscous  terms  (i.e. , -|^  (C  -|^)  and  (£  -|^-))  were  not  split.  Thus 
equation  (6-20)  is  more  consistent  with  the  present  method  than  equa- 
tion (6-19). 

The  same  question  about  splitting  arises  with  the  second  cross- 
derivative  type 


— (5  ^) 

3x  u 3y;i:i 


Applying  the  cell  Integration,  equation  (6-21)  becomes 


(6-21) 


s i «-22> 

Neither  term  is  zero  here  (except  for  the  cell  nearest  corner  point  D 
where  <p  = v and  when  £ « v),  but  whether  they  are  split  or  not,  they 

need  to  be  evaluated.  Since  the  velocity  on  the  wall  is  zero,  the  ques- 
tion arises  (as  in  the  flux  terms)  whether  a simple  form  like  equation 
(6-2)  is  sufficient.  Here  the  answer  is  yes.  Both  equations  (6-2)  and 
(6-13)  were  tried  on  the  back  wall  for  these  viscous  cross  derivatives 
with  no  detectable  difference.  As  an  example  of  the  simple  form,  the 
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first  term  of  equation  (6-22)  becomes 

^x+(-2irL')xf1  “ {T^i+lj+Cij)[4^i+lj+l'hJ,ij+l'h?,i+lj+<l’ij)“0^Ay} 

(6-23) 

The  last  modification  which  was  tried  before  the  actual  cause  of 
the  wiggles  was  discovered,  was  to  make  the  back  wall  a slip  wall.  This 

was  done  by  making  (-^)  nonzero.  The  value  of  v on  the  wall  was  given 

dy  <jj 

by  an  extrapolation  through  the  nearest  two  grid  points.  This  slip  wall 
tended  to  lower  the  v-component  of  velocity  in  the  back  wall  cells  near 
corner  point  C by  up  to  10%  (though  one  point  changed  36%),  and  it  also 
increased  the  pressure  in  these  cells  by  about  5%.  Nearer  the  center- 
line  and  out  into  the  rest  of  the  flow  field  the  effect  was  negligible. 

All  of  these  modifications  had  little  effect  on  the  wiggles,  and 
seemed  to  point  out  the  relative  insensitivity  of  the  solution  to  various 
one-sided  forms  at  some  boundaries.  The  greatest  effect  on  the  solution 
in  the  cells  along  the  back  wall  was  the  difference  between  the  flux 
derivatives  (equations  (6-12)  and  (6-13)).  For  the  upper  wall,  the 
inflow  conditions  were  important  (see  Section  6.1.5). 

For  some  finite  difference  schemes,  a point  like  the  one  nearest 
comer  point  C,  where  the  two  walls  intersect,  can  be  difficult  to 
formulate.  No  real  problem  was  encountered  here  as  the  cell-integration 
method  made  clear  which  terms  needed  to  be  modified.  For  this  cell  (whose 
indices  are  iw.jw),  the  x-  edge  was  considered  normal  to  the  upper  wall 
and  the  y-  edge  normal  to  the  back  wall. 


■ 


* 
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6.1.2  Centerline 

Initially,  all  of  the  finite  difference  forms  used  for  the 
centerline  were  derived  in  a similar  manner  as  those  on  the  wall  bound- 
aries. Oscillations  of  period  2Ay  (better  known  as  wiggles)  appeared 

in  the  solution  for  density  (and  pressure)  normal  to  the  centerline  and 

33 

extended  some  distance  into  the  flow  field.  Kothari  and  Anderson, 
solved  the  Navier-Stokes  equations  for  the  nonreacting  case  of  the  near 
field  of  a chemical  laser.  Their  solutions  of  the  supersonic  flow 
between  the  centerlines  of  two  adjacent  nozzles  had  wiggles  normal  to 
the  centerlines.  They  cited  central  differencing  as  the  cause  of  the 
wiggles.  Central  differencing  was  definitely  not  the  problem  here,  but 
rather  treatment  of  the  normal  derivative.  Consider 

(i^y_ = <1^  " 0 . * * v (6_24) 

For  (t^-)  * (*  , - A ),  a is  zero  for  all  the  flux  terms  because 

3y  ij  Ay  'V-  vy-  V* 

v = 0 on  the  centerline.  But  for  <p  = pe  * p,  an  expression  is  needed. 
Initially,  the  simple  first-order  expression 

(pe)y_  « (Pe) ij  <6"25) 

was  used.  This  resulted  in  the  wiggles.  When  equation  (6-25)  was 
replaced  by  the  second-order  expression 

(pe)y_  “ \ [8(P«)ij  ' (6-26) 


similar  to  equation  (6-16) 
mal  derivative,  (3v/3y)^  , 


, the  solution  became  smooth.  The  nonzero  nor- 
was  based  on  the  second  order  form  from 


..  j 
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equation  (5-9): 


lw»V  " w t9(v)«  ' (VW 


(6-27) 


6.1.3  Downstream  Boundary 

The  conditions  at  the  downstream  boundary,  <p . . , are  not  known 

ii  t j 

before  hand,  hence  some  sort  of  extrapolation  is  normally  employed.  The 
extrapolation  can  either  be  explicit  or  implicit.  Several  schemes  of 
each  were  tried. 

Explicit  forms  are  applied  after  the  interior  points  have  been 
computed  and  are  lagging  the  solution.  The  linear  form 


>ii,j  ’ 2*ii-l,j  " *ii-2,j 


. (6-28) 


where  ii  ■ downstream  column  cell  index,  caused  no  problems  for  the 
time  step  At  * AtCFL>  but  the  quadratic  form 


"ii.j  “ 3*ii-l,j  " 3*ii-2,j  + *ii-3,j 


(6-29) 


resulted  in  diverging  spatial  oscillations. 

When  At  was  increased,  equation  (6-28)  resulted  in  divergence  at 
the  downstream  boundary,  apparently  due  to  lagging.  An  explicit  zero 
gradient  expression 


^ii,  j " 3 [4<tlii-l,j  “ ^ii-2, j ^ 


(6-30) 


was  tried  with  At  » 4 At  . Equation  (6-30)  was  obtained  by  fitting  a 

CfL 


quadratic  through  (ii-2,j),  (ii-l,j),  and  requiring  that  (r*-) 


3x  ii, j 
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Convergent  solutions  were  obtained  but  it  was  later  discovered  that 
wiggles  in  p,  e,  and  p were  created  at  the  downstream  boundary  as  soon 
as  the  shock  began  to  cross  it.  These  wiggles  then  proceeded  to  travel 
upstream  and  increase  in  amplitude  near  the  back  wall,  and  then  decrease 
to  small  amplitude  at  the  downstream  boundary.  This  rather  unusual 
sequence  is  shown  clearly  in  the  3-dimensional  plots  in  Figure  8.  This 
figure  shows  the  pressure  as  a surface  at  six  different  time  steps.  The 
initial  pressure  was  constant  everywhere,  as  shown  in  Figure  8(a).  (The 
wall  is  represented  by  the  zero  values  in  each  plot.)  At  first  (Figure 
8(b))  the  pressure  quickly  dropped  in  the  back  wall  region  while  the 
downstream  pressure  remained  near  its  initial  value.  This  resulted  in 
a recompression  wave  that  travelled  downstream  to  become  the  shock. 

While  the  wave  intersects  the  upper  boundary  (Figures  8(b)  and  8(c)), 
the  pressure  appears  relatively  smooth  except  for  the  inflow  region  (see 
Section  6.1.5).  As  the  wave  moves  to  the  downstream  boundary  (Figure 
8(d)),  the  wiggles  form  normal  to  the  downstream  boundary.  As  previously 
mentioned,  the  wiggles  then  travel  upstream  to  the  back  wall  and  inflow 
regions  and  reduce  in  amplitude  near  the  downstream  boundary.  Notice 
the  inflow  pattern  change  (Figures  8(e)  and  8(f)). 

The  fact  that  the  wiggles  occurred  normal  to  the  boundary  is  typi- 
32 

cal.  Moretti  points  out  that  wiggles  form  in  the  direction  along 
which  there  is  an  inconsistency.  Indeed,  treatment  of  any  boundary 
involves  modification  of  terms  normal  to  it.  Hence  poor  treatment  there 
can  cause  wiggles  which  are  normal  to  the  boundary.  This  also  occurred 


on  the  centerline 
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When  this  was  discovered,  two  implicit  schemes  were  tried. 
The  first  used  the  simple  first-order  algebraic  form 

*ii,j  " ♦ii-l.j 


(6-31) 


This  form  resulted  in  wiggles  also. 

The  second  implicit  scheme  specified  the  value  at  the  x+  edge  of 
a downstream  boundary  cell  by  the  linear  extrapolation  form 


" 2(3*ii,j  ” *11-1, 


(6-32) 


This  also  requires  that 


(it)  . (It) 

v3x'x+  ',3x'x- 


(6-33) 


with  equations  (6-32)  and  (6-33),  the  finite  difference  conservation 
equations  can  be  applied  for  the  dependent  variables  at  (ii,j).  A 
derivative  for  a downstream  cell  becomes,  for  example 


(it)  - _L(4,  _ ) . -L[  ( ~ <t>li~1?l)  _ (^ii.j  J’il-liJ.)  ] 

',3x''ii,  j Ax'vx+  ?x-  Ax1 ' 2 ' v 2 


Ax  (<f>ii,j 


) 


(6-34) 


Note  that  equation  (6-34)  corresponds  to  upwind  differencing,  so  that 
the  outflow  column  of  points  has  upwind  differencing.  This  form 
resulted  in  a smooth  density  and  pressure  variation  and  eliminated  the 
wiggles.  It  should  be  noted  that  a higher  order  extrapolation  could  be 
used  instead  of  equations  (6-31)  and  (6-32).  This  would  result,  however, 
in  the  block  coefficient  matrix  no  longer  being  tridiagonal.  While  this 


• 
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i 


, 


is  not  too  difficult  a task  to  overcome  by  Gauss  elimination,  the  extra 
programming  effort  was  not  done  here  since  the  linear  form  was  satisfac- 
tory. 

6.1.4  Upper  Boundary 

The  conditions  at  the  upper  boundary  are  not  known  beforehand 

either.  However,  the  boundary  does  lie  outside  the  regions  of  major 

viscous  transport  and  is  assumed  to  be  nearly  inviscid.  Thus  a simple 

wave  condition,  as  described  in  Section  5.5,  appeared  to  be  a logical 

24 

choice.  Allen  and  Cheng  obtained  stable  and  realistic  results  using 
this  scheme. 

While  the  simple  wave  condition  as  applied  is  an  explicit  proce- 
dure, it  produced  stable  and  realistic  results  here  even  for  At  = 32  AtCFL 
The  time  lagging  appeared  not  to  hinder  convergence.  Here,  the  explicit 
boundary  condition  is  desirable  since  it  is  easier  to  apply  in  a computer 
program  than  an  implicit  form. 

One  other  scheme  was  tried  for  the  upper  boundary  in  an  investiga- 
tion of  the  spurious  compressions  in  the  inflow  region.  This  alternate 
form  was  similar  to  the  implicit  extrapolation  scheme  which  worked  so 
well  for  the  downstream  boundary  (see  Section  6.1.3).  In  this  case  the 
values  at  the  y+  edge  were  specified  as 


S+"  2 (3*i,jj  ‘ <(>i,jj-l) 

This  results  in  a derivative  for  the  upper  boundary  point  as 

(3y^i, jj  “ 4y^i« jj  _ *i,jj-l) 


(6-35) 


(6-36) 
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This  form  did  not  work  well.  The  pressure,  internal  energy,  and  density 
of  the  upper  boundary  became  very  small  near  the  downstream  end  and  the 
vertical  velocity  component  became  absurdly  large  positive.  The  probable 
cause  can  be  seen  from  the  form  of  equation  (6-36).  When  $ is  positive 
this  expression  is  the  first-order  upwind  difference.  For  <t>  negative, 
equation  (6-36)  is  the  downwind  difference  form,  which  has  been  shown  to 
be  unstable  (see  for  example,  Roache,  p.  69).  Near  the  inflow  region  all 
the  dependent  variables  are  positive.  But  the  corner  expansion  causes 
the  vertical  velocity  component,  v,  to  become  negative  on  the  rest  of  the 
upper  boundary.  This  means  that  equation  (6-35)  is  stable  for  the  short 
distance  that  v is  positive  and  unstable  for  the  remainder  of  the  upper 
boundary.  It  also  did  not  eliminate  the  inflow  compressions. 

Since  the  explicit  simple  wave  procedure  appeared  not  to  cause  any 
problem,  even  for  large  At,  it  was  retained  as  the  method  for  computing 
the  row  of  points  along  the  upper  boundary. 

6.1.5  Inflow  Boundary 

Initially,  as  discussed  in  Section  5.6,  the  inflow  boundary  was 
specified  and  held  fixed.  The  conditions  were  derived  by  first  assuming 
a u-profile  and  solving  for  internal  energy  by  the  Busemann  energy 
integral  of  the  compressible  boundary- layer  equations  and  density  by  the 
equation  of  state,  assuming  3p/3y  ■ 0.  The  normal  velocity  component  was 
then  solved  by  combining  the  continuity  and  the  momentum  boundary-layer 
equations  and  using  the  energy  integral. 

Two  u-profiles  were  used,  corresponding  to  those  used  by  Allen 
24 

and  Cheng.  The  first  was  a linear  u-profile  up  to  the  boundary  layer 
edge  (which  was  one  third  the  base  half  height) . When  Allen  and  Cheng 
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used  this  profile  they  had  v - 0 as  required  by  the  formula  for  v.  The 
second  u-profile  was  the  polynomial 

u(n)  * (14n  - 7n4  + 2n7)  (6-37) 

where  n - y/6  and  6 * .41.  Here,  v was  not  zero. 

When  the  linear  u-profile  was  used  with  v = 0,  a series  of  com- 
pressions in  the  inflow  region  caused- a pressure  rise  of  as  much  as  22% 
of  the  freestream  value  near  the  upper  wall.  By  requiring  the  normal 
pressure  gradient  on  the  upper  wall  to  be  zero  by  using  equation  (6-26) , 
this  was  reduced  to  a rise  of  6%.  The  zero  pressure  gradient  is  unreal- 
istic near  the  expansion  corner,  however,  where  v becomes  50%  of  u. 

34 

Allen  also  encountered  pressure  rises  with  v - 0.  The  next 
modification,  then,  was  an  attempt  to  compute  v at  each  time  step.  Two 
schemes  were  tried.  The  first  used  the  explicit,  backward,  linear 
extrapolation  scheme 


'l,j 


2v 


- v. 


2,j  3,  j 


(6-38) 


This  form  caused  diverging  oscillations  in  v which  fed  into  the  other 
variables.  The  second  scheme  used  the  implicit  form 


v 


1.J 


2,J 


(6-39) 


This  resulted  in  stable  computations  and  reduced  the  pressure  increase 
in  the  compressions  to  2%  of  freestream.  The  zero  pressure  gradient 
was  still  applied  and  v continued  to  Increase  outside  of  the  boundary 
layer  along  the  inflow  boundary.  Because  these  results  were  non-physical, 
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and  because  the  linear  u-profile  has  an  abrupt  change  in  slope  at  the  4 

boundary  layer  edge,  it  was  decided  to  use  equation  (6-37),  which  has 

i continuous  first  and  second  derivatives  at  the  boundary-layer  edge. 

Then  v was  computed  by  the  quadrature  of  equation  (5-51)  and 

held  fixed.  For  the  inflow  cells  outside  the  boundary  layer,  v was 

equal  to  the  value  computed  by  equation  (5-51)  for  n = 1.  When  this 

was  applied  and  the  normal  pressure  gradient  condition  removed,  the 

compressions  were  still  there  and  resulted  in  a maximum  pressure  rise 
34 

of  about  3%.  Allen  pointed  out  that  too  few  grid  points  in  the  bound- 
ary layer  could  result  in  an  inaccurate  representation  of  the  large  gra- 
dients in  density  and  internal  energy  (and  hence  pressure)  near  the 
boundary- layer  edge.  If  this  were  true,  then  a finer  mesh  would  reduce 
the  magnitude  of  the  compressions.  This  was  tried  (see  Section  6.2)  and 
indeed  the  compressions  disappeared.  A single  wave  originating  from  the 
top  of  the  boundary  layer  does  exist  in  the  fine  mesh  solution,  but  it 
does  not  have  the  same  character  as  the  other  irregularities.  It 
appears  to  be  an  expansion  wave,  since  the  pressure  drops  immediately 
downstream  of  it.  As  discussed  ection  6.2,  this  wave  may  be  a 
result  of  the  use  of  3P/3y  * 0 along  the  inflow  boundary,  when,  in  fact, 
the  corner  expansion  is  being  felt  this  far  upstream  of  the  corner. 

6.2  Computational  Results 

In  this  section  the  results  of  the  computational  method  are  given 

24  34 

and  compared  to  those  obtained  by  Allen  and  Cheng  and  Allen.  The 
computations  were  made  for  a laminar,  supersonic  two-dimensional  flow 
past  a corner  (see  Figure  2)  at  ■ 


3 and  Re  ■ 550  with  a boundary- layer 
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height  of  41%  of  the  base  half  height.  The  inflow  profiles  were  those 
given  by  Allen  and  Cheng  and  the  initial  conditions  were  those  described 
in  Section  5.7.  The  mesh  size  was  Ax  ■ 1/6  and  Ay  =*  1/12.  This  was 
coarser  than  the  Allen  and  Cheng  mesh  (Ax  ■ 1/9.75  and  Ay  - 1/19.5)  and 
was  chosen  to  increase  computational  speed  since  fewer  grid  points  were 
needed  to  retain  reasonable  accuracy.  The  check  of  the  coarse  mesh 
accuracy  was  done  by  using  a finer  mesh  for  comparison  (see  Section  6.3). 
The  ratio  Ax/Ay  ■ 2 was  maintained.  The  time  step  size  used  in  this  test 
was  four  times  the  maximum  time  step  for  stability  in  explicit  methods, 
known  as  the  Courant-Freidrichs-Lewy  (CFL)  time  step.  Stable  and  conver- 
gent solutions  were  obtained,  demonstrating  potentially  large  savings  in 
computer  time.  This  is  discussed  in  more  detail  in  Section  6.3. 

Figures  9-20  give  the  results  of  the  computations  using  the  above 
conditions.  Figure  9 shows  the  velocity  vectors,  giving  both  magnitude 
and  direction.  The  arrowheads  are  all  the  same  size  so  that  direction 
can  still  be  seen  for  very  small  velocities.  The  main  features  of  the 
flow  field  can  be  seen,  including  the  upper  wall  boundary  layer  along  AB, 
the  expansion  and  turning  at  corner  point  C,  the  weak  shock  and  turning 
downstream,  the  retarded  flow  near  the  centerline  DE,  and  the  recirculat- 
ing region  near  the  back  wall  CD.  A rear  stagnation  point  on  the  center- 
line  and  the  separation  point  below  the  corner  are  also  evident.  That 

the  flow  separates  below  the  corner  has  been  shown  in  experiments  by 
35  36 

Hama  and  by  Donaldson.  By  linear  Interpolation  to  find  where  v ■ 0 
in  the  column  of  points  nearest  the  wall,  the  flow  appears  to  separate 
20-25%  of  the  base  half  height,  H,  below  the  corner.  The  separation 
point  obtained  by  Allen  and  Cheng  was  roughly  20%  of  H below  the  comer. 
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The  rear  stagnation  point  here  was  2%  closer  to  the  wall  than  the  one 
computed  by  Allen  and  Cheng. 

Some  of  these  features  are  more  clearly  shown  in  the  streamline 
plot  of  Figure  10.  Expansion  and  compression  are  shown  by  the  spreading 
and  converging  of  the  streamlines.  Flow  direction  is  indicated  by 
streamline  slope  (this  also  illustrates  the  inflow  through  AF) . The 
recirculation  region  is  characterized  by  the  closed  loops  near  the 
back  wall  and  the  dividing  streamline  (denoted  by  the  "+"  symbols)  is 
shown  extending  from  the  back  wall  to  the  centerline. 

The  value  of  the  stream  function  was  computed  at  each  point  in 
the  flow  field  by  summing  the  mass  flows  vertically  starting  from  the 
zero  streamline  BCDE.  For  example, 

^ + (Ay/2)[(pu)i;j  + (6-40) 

Maximum  and  minimum  \p  were  found  to  determine  the  range.  The  range  was 
then  divided  into  a number  of  incremental  values.  Between  the  minimum 
value  of  i p and  ii  * 0,  two  equally  spaced  values  were  determined.  Between 
ifi  * 0 and  the  maximum  value  of  , the  increments  were  evenly  spaced  along 
the  inflow  boundary.  With  the  array  of  incremental  values,  the  plotting 
routine  searched  row  by  row  and  column  by  column  to  find  where  these 
incremental  values  occurred  between  grid  points.  Linear  interpolation 
was  used  here  to  determine  the  variation  between  grid  points.  The  loca- 
tion of  the  separation  point  given  by  the  dividing  streamline  is  not 
accurate  for  two  reasons.  The  first  is  that  the  back  wall  region  has 
been  shown  to  be  sensitive  to  the  boundary  condition  treatment , thus  the 
variables  themselves  are  suspect.  Sunning  over  several  cells  adds  all 
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the  errors.  Thus  determination  of  the  separation  point  by  interpola- 
tion between  velocity  vectors  is  thought  to  be  more  accurate.  Second, 
the  location  of  the  dividing  streamline  is  not  given  at  the  wall  but 
half  a cell  away.  Since  the  wall  is  also  a t 5 0 line,  the  angle  of 
approach  is  not  known  exactly. 

The  contour  plots  (Figures  11  through  13)  were  made  in  a similar  way 

as  the  streamline  plot,  except  that  there  were  equal  increments  between 
maximum  and  minimum  values.  The  maximum,  minimum,  and  incremental 
values  of  the  variable  are  pointed  out  above  each  plot  with  the  corre- 
sponding symbols. 

Figure  11  clearly  shows  the  effect  of  the  expansion  and  compression 
on  the  pressure.  As  expected,  the  pressure  drops  rapidly  as  the  flow 
expands  around  the  corner.  Note  that  the  expansion  begins  upstream  of 
the  corner.  The  single  wiggle  in  the  "Y"  line  upstream  of  the  corner 
and  the  wiggles  in  the  " O " line  at  the  downstream  boundary  appear  to 
be  caused  by  the  coarseness  of  the  mesh.  In  the  fine  mesh  solution, 
these  two  wiggles  are  nearly  gone.  A growing  region  of  3p/3y  5 0 near 

% 

the  centerline  indicates  the  region  where  use  of  the  boundary-layer  equa- 
tions would  be  a valid  representation  of  the  viscous  wake.  The  pressure 

distribution  across  the  base  was  relatively  constant  except  for  a sharp 

35 

50 Z drop  near  the  separation  point.  Hama  made  measurements  of  the 
base  pressure  variation  for  ■ 4.54  and  22  x 10^  < Re  < 2.03  x 10^ 
while  studying  the  lip  shock.  His  results  show  the  same  behavior.  He 
found  that  the  lip  shock  strength  increases  with  Mach  number  and 
Reynolds  number,  but  decreases  with  boundary  layer  height.  This  proba- 
bly explains  why  no  lip  shock  was  observed  in  these  calculations.  The 
Reynolds  number  was  only  550,  the  boundary- layer  relatively  thick,  and 
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the  Mach  number  not  particularly  high. 

4 

The  density  contours  in  Figure  12  and  the  corresponding  internal 
energies  in  Figure  13  also  show  the  effect  of  expansion  and  compression. 

Near  the  adiabatic  wall  the  temperature  (e/cv)  is  higher  than  the  free- 
stream  and  the  density  shows  a steep  gradient  in  the  upper  part  of  the 
boundary  layer.  The  linear  nature  of  the  contour  lines  for  pressure 
(Figure  11),  density,  and  internal  energy  near  the  upper  boundary  show 
the  validity  of  the  simple  wave  condition.  The  maximum  density  was 
within  1.5%  of  Allen  and  Cheng,  the  minimum  within  17.0%.  The  maximum 
internal  energy  was  within  1.0%  and  the  minimum  1.4%. 

Figure  14  compares  the  results  for  the  pressure  along  the  center- 

line  with  the  results  of  Allen  and  Cheng.  The  agreement  is  very  close. 

The  Allen  and  Cheng  line  was  taken  directly  from  their  results.  The 

present  results  line  is  not  actually  on  the  centerline  but  Ay/2  above 

it.  That  this  should  be  an  accurate  representation  of  the  pressure  on 

the  centerline  can  be  argued  by  the  fact  that  dp/dy  = 0 on  the  center- 

2 

line.  Hence  the  error  is  of  the  order  of  (Ay)  by  Taylor  series  expan- 
sion. Also,  the  pressure  contours  show  that  there  is  little  variation 
between  the  pressure  on  the  centerline  and  at  points  nearby. 

Further  close  agreement  with  the  pressure  computations  of  Allen 
is  shown  in  pressure  profiles  in  Figure  15.  Here,  the  vertical  pressure 
variations  at  several  x/H  locations  are  plotted.  Once  again  3p/3y 
approximately  zero  near  the  centerline  is  indicated.  The  pressure  drop 
due  to  expansion  is  shown  by  the  decrease  in  pressure  from  the  upper 
part  of  the  plot  down  to  about  y/H  * 1.0.  Then  the  pressure  rises,  more 
sharply  as  X/H  increases,  due  to  the  flow  turning  and  the  formation  of 
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of  the  shock. 

This  behavior  can  be  seen  more  clearly  in  the  three-dimensional 
pressure  plot  of  Figure  16.  This  shows  the  pressure  in  the  entire  flow 
field  as  a surface.  The  height  in  the  three-dimensional  plot  is  the 
pressure.  The  "floor"  of  the  box  around  the  contour  is  the  plane 
defined  by  ABCDEF  of  the  two  dimensional  flow  field  (see  Figure  3). 

The  values  inside  the  wall  are  arbitrary,  of  course,  since  they  are  never 
used  and  have  been  set  equal  to  zero.  The  inflow  region  is  at  the  upper 
left  and  the  outflow  at  the  lower  right  (EF) . The  centerline  DE  is 
also  shown.  The  points  A and  C are  hidden.  The  plot  is  easy  to  gener- 
ate, being  a simple  FORTRAN  CALL  in  the  CALCOMP  plotting  package. 

The  expansion  around  the  corner  ia  very  clear  as  the  contour 
shows  a steep  drop  from  the  inflow  value.  The  downstream  development 
of  the  recompression  shock  as  shown  as  the  growing  difference  in  height 
between  the  pressure  near  the  center  of  the  field  and  the  region  above 
the  centerline.  The  extent  of  this  difference  can  be  seen  on  the  verti- 
cal plane  through  EF. 

As  demonstrated  in  Figure  8,  this  type  of  plot  can  be  extremely 
useful  in  visualization,  not  only  of  flow  features,  but  of  the  propa- 
gation of  disturbances.  In  Figure  16  there  are  a series  of  pressure 
fluctuations  near  the  Inflow  and  at  the  "crest"  of  the  recompression 
wave.  Those  at  the  inflow  are  present  from  the  beginning  of  the  compu- 
tation (see  Figure  8),  while  those  in  the  shock  appear  only  after  the 
shock  has  settled  to  its  final  position.  Both  appear  to  be  caused  by 
the  coarseness  of  the  mesh,  however,  since  these  irregularities  are 
absent  in  the  fine  mesh  solution. 
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The  centerline  Mach  number  distribution  is  given  in  Figure  17 . 

The  extent  of  the  recirculation  region  and  the  downstream  acceleration 
are  indicated.  Slightly  further  downstream  of  the  computation  region 
the  centerline  Mach  number  will  exceed  one.  The  flow  is  already  super- 
sonic a short  distance  above  the  centerline.  No  singularities  asso- 
ciated with  the  Mach  one  (sonic)  condition,  such  as  the  Crocco-Lees 
singularity,  were  encountered. 

The  values  computed  here  were  somewhat  higher  than  Allen's.  At 
x/H  » 4.75  the  difference  was  about  6%.  The  reason  for  the  higher  Mac. 

t 

numbers  on  the  centerline  were  due  to  the  lower  speeds  of  sound,  since 
the  temperature  was  lower.  These  lower  temperatures  are  indicated  by 
the  lower  internal  energies  as  indicated  in  Figure  18.  The  present 
results  for  the  internal  energy  on  the  centerline  were  typically  12-15% 
lower  than  those  computed  by  Allen.  The  reason  for  the  difference  is  not 
known,  although  it  is  presumed  to  be  a result  of  the  difference  in  the 
way  the  centerline  was  treated  as  a boundary.  Allen  uses  the  "reflec- 
tion principle"  rather  than  one-sided  difference  forms.  He  passes  the 
centerline  through  the  center  of  a cell  rather  than  its  y-  edge.  The 
centerline  conditions  are  then  enforced  in  a row  of  cells  below  those  on 
the  centerline  by  setting  the  dependent  variables  there  equal  to  the 
values  in  the  row  of  cells  above  the  centerline  cells.  As  previously 
stated,  the  centerline  conditions  were  enforced  in  the  present  method  by 
using  appropriate  one-sided  forms  at  the  cell  edges  adjacent  to  the  cen- 
terline. 

Away  from  the  centerline,  the  internal  energy  comparison  with 
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Allen  improves.  This  is  seen  in  the  internal  energy  profiles  of  4 

Figure  19.  Thus  the  variation  in  the  centerline  values  appears  to 
have  been  caused  by  the  difference  in  boundary  condition  treatment. 

6.3  Computational  Experiments 

• This  section  discusses  the  numerical  considerations  and  the 

results  of  some  experiments  concerning  grid  size,  time  step  size,  and 
initial  conditions. 

All  computations  were  done  on  the  CDC  6600  of  the  Georgia  Tech 

Cyber-74  computer  system.  The  storage  capacity  of  this  machine  permits  1 

J 

over  5000  computational  grid  points,  though  that  many  were  never  used.  ( 

The  solution  of  the  four  conservation  equations  using  1056  grid  point  in  J 

the  flow  field  ( coarse  mesh)  required  10  CPU  minutes  to  1 CPU  hour,  depend-  , 

ing  on  the  convergence  criteria  (see  Section  6.3-2),  the  size  of  the 
time  step  (Section  6.3-3),  and  on  the  initial  conditions  (Section  6.3-4). 

The  basic  grid  configuration  used  48  points  in  the  x-direction  and 
24  points  in  the  y-directlon.  The  length  of  the  upper  wall  was  4/3  the 
back  wall  height.  This  was  23%  longer  than  used  by  Allen  and  Cheng.  The 
downstream  boundary  was  also  further  from  the  back  wall  than  in  Allen 
and  Cheng,  being  6 2/3  step  heights  away  compared  with  5.13.  For  the 
finer  mesh  experiments,  the  number  of  grid  points  in  each  direction  was 
doubled  with  all  distances  remaining  the  same.  Thus  the  coarse  mesh 
solutions  used  1056  grid  points  in  the  flow  field,  and  the  fine  mesh  had 
4224. 

6.3-1  Fine  Mesh  Comparisons 

To  test  the  accuracy  of  the  coarse  mesh  it  is  desirable  to  make 
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the  same  computations  with  a much  finer  grid.  Since  the  finite-differ- 
ence equations  are  consistent  with  the  differential  equations  (see 
Chapter  III),  greater  accuracy  should  be  expected  as  Ax  and  Ay  get 
smaller.  In  addition,  the  scale  of  variation  of  some  quantities  is  not 
large  compared  with  the  size  of  the  coarse  mesh.  A smaller  mesh  would 
thus  increase  resolution. 

The  mesh  was  made  finer  by  dividing  each  coarse  mesh  cell  into 
four  equal  cells  (see  Figure  20).  Thus  no  fine  mesh  point  coincides  with 
a coarse  mesh  point,  and  to  compare  the  two  solutions,  values  at  the 
four  fine  mesh  points  located  in  a coarse  mesh  cell  were  averaged.  The 
basic  results  of  the  fine  mesh  computations  are  given  in  Appendix  D. 

Two  comparisons  between  fine  and  coarse  grid  computations  were 
made.  The  first  uses  the  form 


$ 

coarse 


^fine 

fine 


(6-41) 


This  form  is  useful  for  (j>  near  or  greater  than  1 and  represents  a per- 
centage change.  For  small  (i.e.,  near  zero),  equation  (6-41)  can 
give  a large  number,  even  for  a small  change  in  0.  In  this  case  a more 
appropriate  comparison  may  be  the  simple  form 


fine 


coarse' 


(6-42) 


Tables  1 and  2 show  the  results  of  using  both  equation  (6-41)  and 
(6-42).  For  equation  (6—41),  is  also  given.  Table  1 gives  the  maximum 
difference  in  the  flow  field  between  the  coarse  and  fine  meshes.  Most 


of  these  maximum  differences  occur  near  the  back  wall. 


AO-AO50  mi  AIR  FORCE  INST  OF  TECH  MR I OHT -PATTERSON  AFB  OHIO  F/0  20/4 

AN  IMPLICIT  FINITE  DIFFERENCE  PROCEDURE  FOR  THE  LAMINAR*  SUPERS— ETC (U) 
197A  R L ROACH 

UNCLASSIFIED  AFIT-CI-70-71  NL 


80 


Table  1.  Coarse-fine  Mesh  Comparisons  - Flow  Field  Maximuas 
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Table  2.  Coarse-fine  Mesh  Comparisons  - Selected  Points 


Point 

U$/$f  I 
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P 

.027 

.167 

.004 

near  C 
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Table  2 gives  the  result  for  some  selected  points  In  the  flow 
field.  Included  are  the  cells  nearest  to  corner  point  C,  corner  point 
D,  corner  point  E,  and  a point  about  4.8  step  heights  downstream  and  In 
the  recompression.  Again,  the  largest  differences  are  near  the  back 
wall.  This  serves  to  point  out  the  sensitivity  of  this  region. 

One  Important  result  of  the  fine  mesh  study  was  to  show  that  the 
coarse  mesh  waa  responsible  for  the  irregularities  near  the  Inflow  and 
those  in  the  shock  near  the  downstream  boundary.  Apparently,  the  finer 
mesh  allowed  greater  resolution.  Figure  21  is  a three-dimensional  plot 
of  the  pressure  using  the  fine  mesh. 

A single  wave  remains  at  the  inflow,  but  its  character  is  differ- 
ent from  those  in  the  coarse  mesh.  This  wave  has  an  angle  of  about  21°, 
relatively  close  to  the  inflow  Mach  angle  (19.5°)  and  appears  to  be 
related  to  the  choice  of  inflow  boundary  conditions.  As  mentioned  pre- 
viously, 3p/3y  was  assumed  to  be  zero  along  the  inflow,  but  it  can  be 
seen  that  the  expansion  corner  is  being  felt  even  this  far  upstream 
(4/3  H)  to  reduce  the  pressure  near  the  wall.  The  wave,  originating  at 
the  top  of  the  boundary  layer,  may  be  a flow  adjustment  because  of  the 
inconsistency.  The  wave  may  also  be  caused  by  the  change  from  the 
boundary  layer  u-proflle  to  the  freestream  u « const,  profile. 

6.3.2  Convergence  Criteria 

As  the  rate  of  change  with  time  of  the  dependent  variables  in 
the  flow  field  decreases,  the  solution  is  said  to  be  converging.  When 
this  rate  is  zero,  the  steady  state  has  been  reached.  It  is  usually 
not  necessary  or  practical  to  continue  the  calculation  until  the  rate 


Figure  21 


Pressure  Surface  - Fine  Mesh 
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is  zero,  however.  Allen,  for  example,  uses  a convergence  criteria  of 


I n+1 
IP 


max 


< e 


(6-43) 


where  2x10  ^ < e < 2x10  ^ . 


In  the  present  work,  the  form  normally  used  was 


i n+1  ni 
|p  - P [max  < 

n 1 

P 


(6-44) 


with  pn  known,  the  range  of  (pDe^)  can  be  compared  directly  with  the  range 
of  Allen's  c: 


2 x 10-6  < pne  < 5.3  xlO-5  (6-45) 

25 

As  in  Allen  and  Cheng,  a check  on  density  was  normally  sufficient  to 
determine  convergence.  Equation  (6-44)  was  applied  also  to  u,  v and  e 
as  an  additional  check. 

The  rate  of  convergence  varied  over  the  whole  flow  field.  The 
most  rapid  convergence  was  in  the  supersonic  regions,  especially  upstream 
of  the  corner.  The  slowest  rate  of  convergence  was  in  the  region  near 
the  back  wall.  It  may  be  possible,  therefore,  to  shorten  the  calculation 
for  cases  where  the  near  wall  region  is  not  a region  of  primary  interest. 
6.3.3  Time  Step  Studies 

As  mentioned  before,  the  chief  advantage  in  using  an  implicit 
procedure  is  the  ability  to  use  larger  time  steps  than  allowable  in 
explicit  methods,  and  hence  reduce  computation  time.  Here,  the  CFL 
stability  condition  in  two  dimensions  is 


4 
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..  . 1 . 1,1.1  xl/2,-1 

AtCFL  1 + ^7(7T+7T)  1 

1 Ax  Ay 


(6-46) 


In  the  present  work,  a time  step  size  of  At  ■ 4At  was  used  while 
developing  the  method.  This  was  chosen  because  it  is  large  enough  to 
be  appreciably  larger  than  At^^  but  small  enough  to  avoid  possible 
problems  with  the  initial  conditions.  Results  for  At  ■ 4 AtCFL  were  given 
in  Section  6.2. 

To  examine  the  ability  of  the  procedure  to  use  larger  time  steps. 

At  ■ 8,  16,  32  and  40  At^^  were  also  attempted.  A time  step  limitation 
was  expected  because  the  finite-difference  equations  are  linearized  with 
respect  to  time  about  the  known  time  level.  This  allows  a simpler  and 
more  rapid  solution  of  the  finite-difference  equations  than  if  they  were 
left  nonlinear.  Still,  the  allowable  At  is  much  larger  than  AtCFL> 
which  more  than  compensates  for  the  larger  computation  time  per  time  step 
relative  to  explicit  methods. 

For  At  ■ 40  AtCFL  the  solution  did  not  converge,  but  diverged  in 
less  than  10  time  steps.  Whether  this  would  happen  with  different  ini- 
tial conditions  has  not  been  determined.  The  rest  of  the  discussion  con- 
cerns the  convergent  solutions  for  At  - 4,  8,  16,  32  AtCFL.  The  differ- 
ence in  the  solutions  at  convergence  for  each  of  these  time  steps  was 
very  slight.  The  difference  in  p between  the  two  solutions  for  At  * 4AtCp^ 
and  At  ■ 16  AtCFL>  for  example  was  typically  around  1Z. 

For  At  ■ 4,  8,  16  At^.^,  the  solution  at  equal  elapsed  times 
appeared  quite  similar,  that  is,  the  solution  in  time  appeared  to  be  rel- 
atively insensitive  to  the  size  of  At.  Table  3 presents  the  solution  at 
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selected  points  in  the  flow  field  for  elapsed  times  of  T - 15.279  and 
T * 76.393  (where  T = H/U^)  for  At  - 4,  8,  and  16  AtCFL.  It  appears, 
then,  that  up  to  16  AtCFL>  for  these  initial  conditions,  doubling  the 
time  step  allows  the  same  point  in  time  to  be  reached  about  twice  as 
fast  and  thus  doubles  the  rate  of  convergence.  This  can  be  seen  in 
Figure  22  which  gives  the  maximum  time  rate  of  change  of  density  of  the 
flow  field  versus  the  number  of  time  steps.  Since  computation  time  is 
directly  related  to  the  number  of  time  steps  taken,  this  is  shown  as  the 
lower  abcissa.  This  figure  clearly  shows  the  substantial  computation 
time  savings  from  the  ability  to  use  larger  time  steps.  The  "humps"  on 
each  curve  are  caused  by  a shift  in  location  of  the  maximum  change  in 
Ap/p  from  a point  near  the  outflow  boundary  to  a point  near  the  back 
wall. 

The  rate  of  convergence  for  At  * 32  At^,^  is  also  given.  Initial 
oscillations  which  damped  out  prevented  the  run  at  this  time-step  size 
from  following  the  same  temporal  path  as  the  other  time-step  sizes.  This 
is  reflected  in  Figure  22  as  the  convergence  was  not  smooth.  When  the 
computation  did  smooth  out,  the  solution  still  converged  faster  than 
16  4tCFL>  though  not  twice  as  fast. 

The  curve  for  At  - ®AtcFL  shows  a sudden  drop  at  t/At  * 150.  As 
indicated,  this  corresponded  to  continuing  the  computation  with  At  ■ 32  A t^,^. 
This  demonstrates  that  convergence  can  be  speeded  even  after  computation 
has  begun  by  increasing  the  time  step  size.  While  no  attempt  was  made 
to  develop  a scheme  for  relating  time  step  size  to  rate  of  convergence, 
such  procedures  can  be  very  useful  to  obtain  very  fast  convergence. 
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The  maximum  time  step  also  seems  to  be  related  to  the  grid  size. 

For  the  fine  mesh  case,  use  of  At  • 32  it—  resulted  in  diverging  oscil- 

Cr  L 

lations  related  to  the  recompression  (see  Figure  23).  However,  the  fine 
mesh  solution  was  stopped  at  t/At  * 85  and  At  was  changed  to  16  At^,^. 

The  oscillations  decreased  and  the  solution  was  stopped  at  t/At  = 300. 
Figure  17  is  the  3-D  pressure  plot.  Convergence  was  greatly  slowed  by 
the  presence  of  the  initial  oscillations.  No  other  time  step  studies 
were  accomplished  with  the  fine  mesh  because  of  increased  run  time  asso- 
ciated with  the  increased  number  of  grid  points. 

As  mentioned  before,  the  computation  time  depends  on  the  size  of 
the  time  step  and  this  was  shown  dramatically  in  Figure  22.  The  number 
of  time  steps  required  for  the  same  degree  of  convergence  was  reduced 
nearly  sixfold  as  At  was  increased  from  four  to  32  At^,^.  In  comparing 
the  computation  time  with  Allen's  method,  it  was  helpful  that  he  had  also 
made  3ome  calculations  on  a CDC-6600  machine.  He  required  about  1.5 
msec  per  time  step  per  grid  point  and  typically  required  between  1000  and 
2000  time  steps  for  convergence.  The  present  procedure  required  7.7 
msec  per  time  step  per  grid  point.  Since  the  present  method  was  able 
to  use  a time  step  30  times  larger,  this  represents  a sixfold  decrease 
in  overall  computation  time  for  these  conditions. 

6.3.4  Initial  Conditions 

To  insure  that  the  given  Initial  conditions  were  not  unique  for  a 
given  set  of  boundary  conditions,  and  to  examine  the  effect  of  Initial 
conditions  on  the  ability  to  take  large  time  steps,  it  is  desirable  to 
examine  several  different  initial  conditions.  Computer  time  and  cost. 


however,  placed  restrictions  on  the  number  of  choices. 


Figure  23.  Pressure  Surface  - Fine  Mesh  with 
At  - 32  AtCFL  aftar  85  time  Steps. 
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Since  the  region  near  the  back  wall  (below  the  corner)  seemed  to 
be  the  most  sensitive  to  other  changes,  a wide  range  of  initial  hori- 
zontal velocity  components  (0  <_  u <_  1.0)  below  the  corner  was  examined. 

The  time  step  for  this  study  was  held  fixed  at  At  - 16  AtrT7T  . 

For  u ■ 0,  the  computation  diverged.  For  0.1  <_  u <_  1.0,  the  compu- 
tations converged.  As  in  the  case  of  different  time  step  sizes,  the 
final  solutions  for  the  different  initial  conditions  were  nearly  iden- 
tical. The  rates  of  convergence,  however,  were  much  different.  Figure 
24  shows  that  for  u • 1.0,  the  rate  of  convergence  was  several  times 
faster  than  when  u * 0.1.  Note  also  that  the  convergence  paths  are  not 
similar.  The  curve  for  u ■ 0.3,  for  example,  has  a rise  at  t/At  ■ 30 
corresponding  to  a shift  of  the  maximum  Ap/p  in  the  flow  field  from  the 
downstream  boundary  to  the  back  wall  region.  The  shift  of  maximum  Ap/p 
to  the  back  wall  for  u - 1.0  occurred  at  t/At  » 62  but  was  not  attended 
by  a rise  in  the  curve.  The  recirculation  region  appeared  to  already  be 
near  its  final  solution  in  this  case. 

6.3.5  Net  Mass  Flux 

The  net  mass  flux  through  the  flow  field  was  computed  at  each  time 
step.  In  the  steady  state,  the  net  mass  flux  should  be  zero.  The  compu- 
tation was  accomplished  by  summing  -puAy  across  the  inflow  boundary  cells 
along  AB,  pvAx  across  the  upper  boundary  cells  along  AF,  and  puAy  across 
the  downstream  boundary  cells  along  EF.  For  coarse  mesh  solutions,  the 
net  inflow  differed  from  the  net  outflow  by  1.6  to  1.8Z.  For  the  fine 
,!  mesh  computations  the  difference  was  1Z. 


Rates  of  Convergence  for  Different  Back  Wall 
Initial  Velocities  At  = 16At„.,,  . 
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CHAPTER  VII 

CONCLUSIONS 

1.  The  results  In  Section  6.2  show  that  the  method  successfully 
computed  all  the  main  features  of  the  flow,  Including  the  corner  expan- 
sion, the  recompression  shock,  the  recirculation  region,  the  viscous 
wake  near  the  centerline,  and  the  simple  wave  nature  of  the  flow  near 
the  upper  boundary. 

2.  Great  care  must  be  taken  in  the  formulation  of  the  boundary 
conditions  to  achieve  physically  realistic  results,  convergence,  and  to 
avoid  wiggles  in  the  steady-state  solution.  First-order  forms  for  pres- 
sure and  3v/3y  on  the  centerline,  for  example,  gave  a converged  solution 
but  had  y-direction  wiggles  in  the  steady  state.  Second-order,  one- 
sided forms  removed  the  wiggles.  For  the  outflow  boundary,  all  explicit 
extrapolation  schemes  caused  divergence  for  At  > AtCFL.  Zero-gradient 
forms,  both  explicit  and  implicit,  gave  x-direction  wiggles  for  the 
converged,  steady  state  solution.  A new  implicit,  linear  extrapolation 
scheme  which  uses  the  finite  difference  equations  was  developed  for  the 
downstream  boundary  and  gave  a smooth,  converged  solution. 

3.  No  artificial  viscosity  was  required  for  stability  and  con- 
vergence. Briley  and  McDonald,  however,  required  additional  explicit 
artificial  viscosity  in  their  subsonic  duct  flow  solutions.  The  reasons 
for  this  difference  are  not  known.  It  may  be  speculated,  however,  that 
the  difference  arises  from  the  present  use  of  the  conservative  form  of 
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the  conservation  equations  (Briley  and  McDonald  use  the  nonconservation 
form  of  the  energy  equation),  the  cell  Integration  technique  for  gener- 
ating finite-difference  equations,  and  the  corresponding  careful  treat- 
ment of  the  boundary  conditions. 

4.  Three-dimensional  contour  plots  were  an  important  diagnostic 
tool.  It  was  not  discovered  that  the  x-dlrection  wiggles  were  caused 
by  the  treatment  of  the  downstream  boundary  until  the  3-D  plots  were 
made.  The  plots  clearly  revealed  that  as  the  recompression  wave  crossed 
the  downstream  boundary,  the  wiggles  formed  and  propagated  upstream  to 
the  back  wall  and  inflow  regions.  Up  to  then,  the  wiggles  were  thought 
to  have  been  caused  by  ill-treatment  of  the  back  wall  boundary  conditions 
or  by  the  cell  Reynolds  numbers  greater  than  two.  The  use  of  upwind 
differencing,  artificial  viscosity,  or  a much  smaller  Ax  were  avoided 

as  they  were  considered  to  be  undesirable  remedies. 

5.  The  results  for  the  contour  plots  showed  qualitative  agreement 
with  Allen  and  Cheng  and  with  Kronzon,  et  al. , and  close  quantitative 
agreement  where  comparisons  were  possible.  The  centerline  pressure  plot 
showed  very  close  quantitative  agreement  with  Allen  and  Cheng.  As  a 
further  check  on  accuracy,  overall  mass  balances  were  computed  at  each 
time  step.  For  the  coarse  mesh  solutions,  the  net  mass  Inflow  rate 
differed  from  the  net  mass  outflow  rate  by  1.6  to  1.8Z.  For  the  fine 
mesh  solutions  the  difference  was  about  1.0Z. 

6.  For  one  set  of  initial  conditions  a time  step  limitation  was 
established  at  about  32  At^^.  The  limitation  was  expected  because  the 
equations  are  linearized  with  time,  even  though  the  method  is  implicit. 
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Compared  with  Allen's  procedure,  this  method  had  a computation  time  per 
time  step  per  grla  point  approximately  five  times  longer,  but  could  take 
five  steps  over  30  times  larger.  This  represents  a six-fold  decrease 
in  computation  time.  In  addition,  the  ability  to  change  (increase)  the 
size  of  the  time  step  during  computation  to  reduce  the  computation  time 
was  demonstrated.  Thus  a time  step  strategy  might  be  successful  wherein 
smaller  At's  were  used  at  the  beginning,  followed  by  increasing  At  as 
the  steady  state  is  approached.  This  would  be  appropriate,  for  example, 
when  the  assumed  initial  conditions  were  very  far  from  the  steady  state 
solution.  Hence,  the  method  appears  to  offer  significant  time  savings. 

7.  The  steady  state  solution  was  quite  insensitive  to  the  choice 
of  initial  conditions,  but  the  time  to  convergence  appeared  to  be  highly 
dependent  on  them.  A range  of  initial  horizontal  velocities  were  applied 
in  the  region  below  the  expansion  corner,  while  the  boundary  layer  on 
the  upper  wall  ahead  of  the  corner  and  the  freestream  conditions  for  the 
rest  of  the  flow  were  the  same  for  these  tests.  It  was  shown  that  an 
initial  u ■ 0 below  the  corner  led  to  divergence  for  At  • 16  At^,^.  As 

u was  increased  from  10%  to  100%  of  the  freestream  value,  increasingly 
faster  times  to  convergence  were  realized.  In  addition,  convergence  was 
shown  for  a significant  range  of  initial  backwall  u. 

8.  Accuracy  of  the  coarse  mesh  results  was  shown  by  comparisons 
with  the  fine  mesh  solution.  Both  solutions  were  in  close  agreement. 

Small,  irregular  disturbances  in  the  inflow  region  and  In  the  shock  near 
the  outflow  boundary  occurred  in  the  coarse  mesh  solutions.  These  are 
attributable  to  the  lack  of  the  resolution  in  the  coarse  mesh  in  the  inflow 
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boundary  layer  and  In  the  shock  at  the  outflow,  as  they  disappeared  In 
the  fine  mesh  solution. 

9.  These  numerical  results  served  to  demonstrate  that  this  pro- 
cedure produced  stable,  convergent,  and  accurate  solutions,  without  the 
use  of  artificial  viscosity,  when  applied  to  this  complex  problem.  To 
the  author's  knowledge,  no  other  implicit  scheme  has  been  successfully 
applied  to  the  multidimensional,  non-linear  Navier-Stokes  equations  for 
the  supersonic  base  flow  problem. 
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APPENDIX  A 

THE  LINEARIZED  CONSERVATION  EQUATIONS 

Conservation  of  Mass 

Also  given  as  equation  (3-1) 


n+1 


- P 


n 


At 


» 3 , n+1  n.  n n+1  n n,  3 , n n+1  n+1  n n n, 

-)  ■ - 3^  (P  u+pu  - P u ] - -^-[p  v +p  v - p v ] (A-l) 


Conservation  of  x-Momentum 

Application  of  the  linearization  to  equation  (2-12)  gives: 


n+1  n , n n+1  _ n n 

n u 

At 


,P “ "V*  + p “u“ ' “ - 2p "u*\  _ 3 r.n+1  n . „_n  n n+1  „_n  n , 1 ,_n+l_n 

( Ji ) * - g^[p  u + 2p  u u - 2p  u + — 5-(p  e 


2 
*5 


n n+1  n n.  1 


+ pe  - P e ) - — cn+l,  3 r nf  1 n n n n+1  n n n n+1 
7 Re  S ] - T-[p  uv+pu  v+puv 
XX  a y 


- n n n 1 en+l, 
2p  u v - — S ] 
Re  xy 


(A-2) 


c A 3u  2 3v 
xx  “ 3 3x  " 3 3y 


s -|a  + |i 

xy  3y  3x 


Conservation  of  -/-Momentum 

Application  of  the  linearization  to  equation  (2-13)  gives: 
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n+1  n n n+1  „ n n 

,p  v + p v - 2p  v , 3r  n+1  n n . n n+1  n . n n n+1  _ n n n 

\ r A 4.  r ' ) ■ • In  II  W n w + ft  II  w — 7rt  (t  ir 


At 


9x 


r lx  * x u u . ii  ir»  x ii  , n n uTi  _ n u n 
Ip  UV+pU  V + p U V - 2p  U V 


2 2 

- ~ s^1]  - f-[p"V  + 2»VvW1  - 2p"vn  + iij  (p^V-f  A"*1 

,M1 


n n 1 _n+l, 
~ p 6 ^ ~ Re  Syy 


CA-3) 


s - 4 3v  2 3u 
yy  3 3y  " 3 3x 


Conservation  of  Energy 

Equation  (2-14)  becomes: 


2 2 

_1_  r n+1  n n n+1  . n n n+1  K . n n . , nr  n+1  n,  n+1  n 

l P e+pe  - 2p  e + p (u  +v  ) + P u u+v  v 

2 / 2 . 2. ,,  _ 3 «■  ,„n+l  n n , n n+1  n , n n n+1  „ n n n. 

-yCu+v).]}*-  — ty(p  ue+pu  e+pue  -2pue) 

2 2 2 2 
, K , n+1  n-nn^.-n  n+1  n n n+1  n . „ n n n n+1 
+ Ip  u(u  +v  )+3pu  u + p u v + 2p  u v v 

2 2 

,n„n,n  , n , . V 3 n+1,  3 , , n+1  n n , n n+1  n 

3P  u (u  + v )]"p^3^e  } - gp-  {y(  p v e +p  v e 

2 2 2 

, n n n+1  n n n,  K r n+1  n,  n . n , , n n+1  n , _ n n n n+1 

+ pve  - 2p  v e ) + j [p  v(u  +v  )+pu  v +2pvuu 

+ SpVV1  - 3 W2  + vn2)]  - en+1)  + K[|j  (unr^x 


■ n n , , 3 , n n t n , , 
+ V T ) + T—  (U  T + VT  ) 1 

xy  3y  xy  yy  J 


t , Tw.  and  T„„  ar«  given  in  Chapter  II. 
**  */  yy 


(A- 4) 
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APPENDIX  B 

FINITE  DIFFERENCE  FORMS  OF  THE  CONSERVATION  EQUATIONS 

Conservation  of  Mass 

Also  given  as  equation  (3-6) 

3p/3t  - - 5x(pu)  - <5y(pv)  (B-l)  \ 


Conservation  of  x-Momentum 

3pu/3t  - - 6x[pu2  + (-iy)pe  - x ] - 5y[puv  - x ] (B-2) 

**  7 

Conservation  of  y-Momentum 

3pv/3t  - - 5x[puv  - x ] - 6y[pv2  + (-^y)pe  - x ] (B-3) 

xy  ym^  77 

Conservation  of  Energy 

[pe  + j (u2  + v2)]  - - 6x{pu[ye  + j (u2  + v2)]  + (p^)  9X 

- (^)(uTxx  + VTyy)}  ‘ 6y{pv[ye  +f  <u2  + v2)l  + <pfe>9y 

- (£)(uTxy  + VTyy>}  (B_4) 

and  «e  given  in  Chapter  II. 

Jtx  yy  x y 
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APPENDIX  C 


ALTERNATING  DIRECTION  IMPLICIT  FORMS  OF 
THE  CONSERVATION  EQUATIONS 


Conservation  of  Mass 


(?  ~ p-)  ■ - 5x[p*u  + pu*- pu]  - 5y[pv] 

At 

** 

n - pv  _ * * **  i 

(* r — ) * - <5x[p  u+pu  - pu]  - 5y[p  v+  v - pv] 

At 


Subtracting  (C-l)  from  (C-2)  a simplified  system  is: 


(C-l) 

(C-2) 


^~At  " " _<sx(p  u + U*  - pu]  - 6y[pv] 


0 


**  * 

£ -P 


**  ** 


At 


■)  “ -<5y[p  v + pv  - 2pv] 


(C-3) 


(C-4) 


Here  a quantity  with  no  superscript  is  considered  as  being  an  n-level 
quantity.  All  the  remaining  equations  are  written  in  the  simplified 
form. 

Conservation  of  x-Momentum 

( p— U ■ +- * -5x[p  u^+2fxiu  - 2pu^  + (— ^y)  ( p e + pe  - pe) 

At  vi/ 


- ^ (3  <5x  u*  - 6y  v)]  - 6y[puv  - Txy] 


(C-5) 


(p*\+  pu**-_P*u-  _ _ 5y[  p**uv  + pu**v  + puv**  - 3puv  - ^ 6y(u**-u)] 

(C-6) 
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Conservation  of  y-Momentum 


* * „ 

,p  v-t-  pv  - 2pv. 

^ At  ' 


* * * 1 * 

- 6x[p  uv  + pu  v + puv  - 2puv  - — (6yu  + 5x  v ) 


- <5y[pv2  + (-^)pe  - x ] 
ymx  77 


**  **  * 


(P  v+ev 


At 


**  l ** 

+ pe  - 2pe)  - — 6y(v  - v)] 


Conservation  of  Energy 


1*  * * K 2 2 **322 

— {p  e+  pe  - 2pe+p  -^(u  +v  ) + pK[u  u + v v - ■=■  (u  +v  )]}  = 
At  2 L 


* * * K * 2 2 2*  *2 

-6x{y(p  ue+pu  e+pue  - 2pue)  + [p  u(u  +v  ) + 3pu  u + pu  v 

+ 2puw  - 3pu(uZ  + v2)]  - 5xe*}  - <5y{pv[ye  + y(u2  + v2)] 

- - 6y  e}  + K[|— (ut  + vt  ) + -I— (ux  + vx  )]  } 

PrRe  } 3x  xx  xy'  3y  xy  yy  J 


{(p**- p*)[e  + -|(u2  + v2)]  + pK[ (u**  - u )u  + (v  -v  )v] 

+ p (e  * - e*) } ■ - <Sy{  (p  * - p)v[ye  + y (u2  + v2)]  + p(v  - v)[ye 
+ |(u2  + v2)]  + pK(u  - u)uv  - 5y(e  - e)} 


x and  x are  given  in  Chapter  It. 
xy  yy  6 


(C— 7 ) 


1 

(C-8) 


(C-9) 


(C-10) 





APPENDIX  D 


COEFFICIENT  MATRICES 


The  coefficient  matrices  for  the  x-sweep  equations  are: 
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A 

At 

0 

0 
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At 
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V 
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e +2-(u  +v  ) 
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1 r n2  en  . 
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*1 
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2 Ax 
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u (ye  -hr(u  +v  )) 
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The  coefficient  matrices  for  the  y-sweep  equations  are: 


n n 
-u  v 


n2  n 


(Ye11  + ^(u^v2))  | -pnKunvn 


n.  n K.  2,  - 2.x 
-p  (ye  +y(u  +3v  )) 
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The  right  hand  side  column  vector  for  the  x-sweep  equations  is: 


aT  + ^(P^11)  - <5y(pnvn) 


- n n 2 n n 

2o  u , n n . o e 


^ r..ro  n n , p e 2 „ n,  .rnnn  1 n..  nx , 

+ (5x[2p  u +£_ 2 — 35^  6yv  J ~ 5y[p  u v (:yu  + ix v )J 


2°3^“+  6x[2pnunvn  + j^  5yun] -6y[pnvn  + 5-| — j^(py  v° -■|<5xun)] 

*1 


2 2 

n/o  n . 3K  , n ^ n . * „ 

p (2e  +~  (u  +v  ))  2 2 

. ££ + 6x[p  u (2ye  +-y-(u  +v  ))]  - 

i 2 2 

\-5y[pnvn(Yen  + |(un  +vn  ))"p^  <Syen]  + KUxCut^  + vt^) 


+ <5y(ur  + vt  ) } 
xy  yy 


The  two  right  hand  side  matrices  for  the  y-sweep  equations  are: 


* n . n * 
> u + p u 


* n . n * 
) v + p v 


",  n , i,  2. . . a,  " n^*  n.  , n * 
p (e  +pu  +v  ))  + pl(u  u +v  v )+p  e 
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APPENDIX  E 

FINE  MESH  RESULTS 

The  following  figures  are  the  results  for  the  fine  mesh  compu- 
tations. There  were  4224  grid  points  in  the  flow  field.  The  ratio  of 
Ax/Ay  was  2.0  with  Ax  ■ 1/12.  The  inflow  Mach  number  was  3.0  and  the 
Reynolds  number  was  550.  The  time  step  size  used  to  reach  convergence 
was  16  AtcpL. 
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Figure  E-l.  Velocity  Vectors  - Fine  Mesh. 


Streamlines  - Fine  Mesh 


Density  Contours  - Fine  Mesh 


112 


I 


I 

,5 


CD 
vr  ® 

X c*> 

CM 


. LO 

>-  in 

CM 


cn 
IX  m 

m 

CM 


in 
4-  o 

CM 


. cn 
•0  co 


cn 
X in 


+ 


cn 

cn 


O 


to 

CO 

CO 


o 


CO 

_ •«* 
□ CO 

o 


■fc 


- 


« 


Figure  E-5.  Internal  Energy  Contours  - Fine  Mesh. 
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